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A Study of Flow in Stratified Reservoirs by Use 


Of the Thermal Analogy 


CHARLES H. PICKERING* 
N. T. COTMAN* 

JUNIOR MEMBER AIME 
PAUL B. CRAWFORD 
MEMBER AIME 


ABSTRACT 


A heat-conduction model has been developed to 
study the flow of fluids in a stratified oil reservoir 
which is being subjected to unsteady-state deple- 
tion. To simulate stratification, plates of different 
metals were joined together so that cross flow could 
occur between strata. 

A description of the model, plate construction and 
necessary instruments is included. The data indi- 
cate that stratified oil reservoirs experiencing cross 
flow may be treated as a single, homogeneous, pro- 
ducing sand having properties intermediate between 
those of the layers making up the system. A tech- 
nique for averaging core-analysis data to arrive at 
the proper average value of rock permeability and 
porosity is presented. 


INTRODUCTION 


Although oil has been produced from reservoirs 
for 100 years and reservoir injection processes have 
been undertaken for 40 years, the industry is still 
inadequately informed on the problem of reservoir 
inhomogeneities. 

Stratification is perhaps the best-known type of 
reservoir heterogeneity. Uren! discussed the possi- 
ble significance of this nonuniformity in his 1927 
paper on the theoretical aspects of water flooding. 
Since this time the problem of vertical permeability 
variations has been the subject of much discussion 
and numerous papers. 

Levorsen? simply describes the origin of layered 
beds in his Geology of Petroleum. Some beds are 
field-wide and may be traced from well to well, while 
others apparently pinch-out and become discontinu- 
ous. In addition, the physical characteristics of a 
layer may vary from place to place. For these rea- 
sons, a stratified system is not easily defined. 


In calculations involving oil reservoirs in which 
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permeabilities vary, methods presently used con- 
tain many simplifying assumptions. The reservoir 
frequently is considered as consisting of permeable 
layers, each of uniform vertical and lateral permea- 
bility and each behaving as if separated from adja- 
cent layers by thin impermeable layers. Other as- 
sumptions frequently made in reservoir calculations 
are steady-state flow conditions, unit mobility ratio, 
linear geometry, average or homogeneous porosity, 
and negligible capillary-pressure and gravity effects. 

Muskat? has developed solutions to the stratifi- 
cation problem for mobility ratios other than unity 
for both linear and exponential permeability distri- 
butions. He assumes lateral uniformity and continuity 
of all productive strata and does not consider poros- 
ity. 

This study will take into account transient effects, 
cross flow and varying porosity, as well as the 
variations in permeability, by means of a heat 
analogy. 

Muskat> has shown the similarity of the equations 
which describe the flow of fluids in a porous medium 
to the flow of heat through a solid. Landrum, et al,® 
made an analogy between the flow of fluid in an 
unsymmetrical reservoir and the conduction of heat 
in a similarly shaped metal plate. This analogy is 
summarized briefly as follows. If fluid is produced 
from an oil reservoir above the bubble point, the 
resulting pressure distribution will be identical to 
the temperature which would result in a metal or 
heat-conducting plate if heat were similarly removed. 
This means that a study of heat flow in stratified 
thermal models will behave like fluid flow in strati- 
fied oil reservoirs. For mathematical correspond- 
encies, see Ref. 6. 

The fluid mobility &/y is proportional to the thermal 
conductivity, and the porosity-times-compressibility 
product is proportional to the density times specific 
heat of the thermal plate. 


EXPERIMENTAL EQUIPMENT AND PROCEDURE 


The stratified reservoir model used in this study 
is illustrated schematically in Fig. 1. A Leeds and 
Northrup Speedomax Type G temperature recorder 
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was used in conjunction with sixteen 36-gauge 
copper-constantan thermocouples to measure and 
record temperature distribution in both uniform and 
stratified thermal models. This instrument records 
a temperature every two seconds. 

The thermocouples were spaced laterally along 
both sides of the plate in such a manner as to give 
the temperature distribution throughout the plate. 
The change in temperature with time at various 
points in a linear copper-Wood’s metal plate is shown 
in Fig. 2, The parameter X is the ratio of the dis- 
tance of a point measured from the end away from 
the source, to the length of the plate. Thus, the 
curve for X =.9 shows the change in temperature at 
a point close to the source, while that for X = 0.0 
designates temperature data taken at the exterior 
boundary of the plate. 

The linear plates were rectangular, 1-ft long, 4-in. 
wide and of varying thicknesses. The metals used 
were chosen according to their thermal properties 
so that studies could be made of the effects of large 
changes in strata permeability and porosity. 

In the construction of the copper-lead plate, a 
plaster mold was used. The copper plate was placed 
into the mold and heated. One 12- x 4-in. side was 
covered with a thin layer of 50:50 lead-tin solder. 
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FIG. 2— THE CHANGE IN TEMPERATURE WITH TIME 


IN A STRATIFIED, LINEAR, COPPER-WOOD’S METAL 
PLATE. 
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Molten lead then was poured onto this surface and 
allowed to cool. A good bond was achieved in this 
manner. Wood’s metal was bonded to brass and copper 
in much the same way, except that the brass and cop- 
per surfaces were prepared by use of soldering paste 
and Wood’s metal rather than solder. The aluminum. 
Wood’s metal plate was more difficult to construct be- 
cause the Wood’s metal would not adhere to the alumi- 
num. The problem was solved by drilling four 1/8-in, 
holes through the aluminum plate in a rectangular 
pattern. The molten Wood’s metal was then poured 
onto the plate, filling the holes. Upon cooling the 
plates were found to be closely bound together. One 
three-layered plate was constructed by bonding 
Wood’s metal to the copper side of a previously 
built copper-lead plate. Finally, a layer of Wood's 
metal was added to the lead side to make a four- 
layered plate. Other three-layered plates included 
a layer of rubber sandwiched between two metals 
held together by machine screws. 

The plates were mounted with one end in contact 
with a mixture of ice and water in an ice bath, as 
shown in Fig. 1. The temperature of this end was 
maintained constant by directing a vigorous flow 
over the plate edge by means of a stirrer. This 
technique simulates reservoir depletion with a con- 
stant pressure maintained at the producing sand 
face. 

The plate was surrounded by a Lucite cylinder, 
which in turn was wrapped with aluminum foil to 
retard heat gain through radiation. A vacuum pump 
was used to reduce convection effects. 

After the recorder was activated and all points on 
the model had attained the same temperature, a run 
was commenced by suddenly introducing ice and 
water into the bath, which corresponds to an instan- 
taneous drop in reservoir pressure at one end of a 
linear reservoir. Since the recorder registered a 
temperature every two seconds, a temperature at 
each of the 16 thermocouples was recorded once 
every 32 seconds. 

Runs were made on six plates consisting of one 
metal only (simulating uniform homogeneous reser- 
voirs), on four bimetallic plates, on one three-layered 
plate and on one four-layered plate. For the multi- 
layered system, the plates were in direct contact, 
representing a reservoir with no barriers to flow 
between strata. In addition, one plate consisting of 
layers of aluminum, rubber and steel, and two con- 
sisting of layers of brass, rubber and steel were 
investigated. These latter studies correspond to 
depletion when the permeabilities differ by about 
500:1 (see Table 1). 

The temperature-distribution data obtained from 
the runs were plotted as dimensionless groups in 
the manner of Carslaw and Jaeger.’ The equations 
presented by these authors considered heat losses 
to be nonexistent and were not applicable in this 
study. Corrections were made for heat gains through 
radiation and for resistence to heat flow caused by 
‘*skin’’ effects which would apply to that area of 
the plate in contact with the ice water.8 
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The temperature distribution at some real time 
would plot as a curve which would correspond to a 
dimensionless time. Knowing the real time and 
length of the plate, effective diffusivity was calcu- 
lated easily. From correlations involving this param- 
eter, the behavior of the flow of heat in the stratified 
system used is elucidated. By analogy, the flow of 
fluids in layered systems can then be inferred. 

To better evaluate the multilayered systems, runs 
were made on the individual metal plates. These 
data are compared to calculated diffusivities in Fig. 
3 and Table 1. With the exception of steel, the ex- 
perimental data were slightly higher than the cal- 
culated diffusivities. Since the handbook diffusivity 
data are also experimental, it is felt that the meas- 
urements made in this study are valid. 

The results taken from four bimetallic plates are 
shown in Fig. 4. All four show that the diffusivities 
decrease with time during the early production 
period. This may be explained on the assumption 
that, initially, the larger portion of the heat removed 
from the model comes from the more-conductive plate. 
The less-conductive plate is also losing heat, but 
at a lower rate. The combined heat flow would then 
indicate a higher-than-average diffusivity. As the 
temperature wave moves toward the outer boundary, 


the diffusivity decreases until a_ substantially 
steady-state condition is achieved. The diffusivity 


of a system depends on the diffusivities of its 
single-layer components, as shown by the positions 
of the copper-lead and copper-Wood’s metal lines. 
The slopes of the two lines differ but little, while 
the effective diffusivity is lowered considerably due 
to the low conductivity of Wood’s metal compared 
to that of lead. 

The thickness of the plates also causes a shift- 
ing of the curves as seen in Fig. 5. The copper-lead 
plate was converted to a three-layered system by 
bonding a layer of Wood’s metal on the lead side of 
the plate. The results show a downward shifting of 
the curve, as well as a reduction in its slope. 
Another layer of Wood’s metal, attached to the cop- 
per, further lowered and flattened the curve. This 
curve could be considered as having negligible 
slope, indicating a constant diffusivity after a fairly 
short period of time. 

Having the measured diffusivities of several 
layered systems, it was desired to find a satisfac- 
tory method for calculating the diffusivities of the 
same systems. In doing so, two methods were used, 
both of which are weighted diffusivities as follows. 
n 
oF kn on 


—_ n 


joey 
= (C,, Pn by, 
n=1 
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TABLE 1 — PHYSICAL PROPERTIES OF SOME METALS, ALLOYS AND RUBBER 


Btu/(hr)(sq ft)(° F/ft) 


Conductivity, 
Diffusivity, sq ft/hr 
































Specific* Density,* 

Material Heat lb/ft 
Aluminum 0.22 168.5 117 
Brass 0.094 534.4 56 
Copper 0.0921 555 224 
Lead 0.0316 687 20 
Stainless Steel 0.107 488 16 
Wood's Metal 0.0352 659.23 74" 
Rubber 0.109 

*From Ref. 11. 

**From Refs. 9 and 10. 
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FIG. 3—COMPARISON OF CALCULATED AND EXPER- 
IMENTAL DIFFUSIVITIES FOR VARIOUS METAL 
PLATES. 
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Experimental Calculated Experimental 
120 3.15 3.23 
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FIG. 4 — EFFECTIVE THERMAL DIFFUSIVITIES OF 
BIMETALLIC PLATES. 
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where k, = conductivity of the nth layer or k/p 
from core analyses, 

b, = thickness of the nth layer or thickness 
having permeability k from core 
analysis, 

Cm = specific heat, 

p = density, 
pC,, = density-specific heat product equiva- 


lent to product of rock porosity and 
fluid compressibility, 
Nn = diffusivity of the nth layer in thermal 
or fluid-flow units, and 
7, 7 = effective diffusivity in either system. 

Based on analytical studies, Eq. 1 should be 
expected to provide a fairly good correlation of the 
experimental data. 

The values calculated for the plates used are 
shown in Table 2 and are compared to the diffusiv- 
ities obtained experimentally. The values used for 
measured diffusivities were obtained by computing 





TABLE 2 — COMPARISON OF CALCULATED AND 
EXPERIMENTAL DIFFUSIVITIES 


Effective Diffusivity @, (sq ft/hr) 

























n n 
x ah 
a ee 
. < E i tal 
xperimenta 
Plate 2 (Cm Pn Pn x An (Average) 
Aluminum- 

Wood's Metal 2.17 1.83 2.24 
Brass-Wood's Metal 0.89 0.73 1.07 
Copper-Lead 3.54 2.88 3.51 
Copper- 

Wood's Metal 2.98 2.20 2.97 
Copper-Lead- 

Wood’s Metal 3.17 2.45 3.20 
Wood's Metal-Copper- 

Lead-Wood’s Metal 3.11 2.37 3.08 
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FIG. 5 — THE CHANGE IN DIFFUSIVITY 
WITH TIME OF MULTILAYERED, LINEAR 
METAL PLATES. 





the area under the curves, such as Fig. 5, and 
dividing by the total time elapsed. These values 
are in close agreement with diffusivities calculated 
by Eq. 1, with the exception of the brass-Wood’s 
metal plate. 

Further confirmation of the equation was obtained 
by studying a stratified model composed of 0.104 in, 
of rubber between two 1/8-in. copper plates. The 
thickness of the rubber layer was increased to 0.208, 
0.312, 0.392 and 0.56 in. for successive runs. This 
resulted in ratios of total rubber thickness to total 
copper thickness of 0.416, 0.832, 1.248, 1.568 and 
2.240, respectively. Experimentally determined val- 
ues of the rubber used in this experiment indicated 
a heat capacity of 0.6 that of copper and a diffusivity 
of about 0.002 that of copper. Comparison of the 
calculated and the effective diffusivity values of 
the stratified models are shown in Table 3. The 
effective and calculated values agree quite well. 

These results are plotted as ratios of effective 
to calculated diffusivity vs dimensionless time in 
Figs. 6 and 7. These plots then may be used for 
correcting a calculated diffusivity of some system 
to its effective diffusivity at any time within the 
range of dimensionless times studied. Experimen- 
tally, the two-layered systems showed a greater 
change in diffusivity with time than the multilayered 
plates studied, as indicated by the slopes of the 
lines. Mirror-image conditions should be kept in 
mind when building and operating the model. 


DISCUSSION 


The application of the results of Eq. 1 may best 
be illustrated by scaling the models to a hypotheti- 
cal oil reservoir operated at pressures above the 
bubble point, and computing the cumulative fluid 
influx by the Hurst-van Everdingen equation.* 


ca 
7 pL 


+[P(t;)) — P(t2)] Vt, — ty te. coe t 


Q(tm) = 2A (Po — PCti)) Vim 


(3) 


where Q = cumulatéve fluid influx, 


A = cross-sectional area perpendicular to 
the direction of flow, 

C = compressibility of fluid, 

od = weighted porosity, 





TABLE 3 — COMPARISON OF EXPERIMENTALLY 
OBTAINED DIFFUSIVITIES WITH THOSE CALCULATED 
BY PROPOSED EQUATION 








Effective Effective Diffusivity from 
Rubber-Copper Diffusivity Experimental Data 
Thickness Calculated by  — *s = 

Ratio ProposedEq  1%=0.195 #=0.391 t= 0.782 

416 3.52 3.26 3.19 3.14 

-832 2.94 2.92 3.02 2.95 

1.248 2.52 2.36 2.69 2.52 

1,568 2.27 2.24 2.24 2.30 

2.240 1.88 ~ 1.79 1,96 
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FIG. 6 — CHANGE IN DIFFUSIVITY WITH TIME FOR 
BIMETALLIC PLATES. 
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FIG.7—THE CHANGE IN DIFFUSIVITY WITH TIME FOR 
MULTILAYERED PLATES. 


initial reservoir pressure, 
= pressure, 

weighted permeability, 
viscosity, 


I AIO 
" 


total period of time under consideration, 

and 

t = time. 
Because this equation is applicable to unsteady- 
state flow in infinite linear systems, a time will be 
chosen such that the pressure wave has not yet 
reached the outer boundary of the reservoir. The 
equation must be simplified for these calculations 
due to the manner in which the data were obtained. 
The sudden application of ice water to the edge of 
the plate is analogous to a sudden pressure drop at 
the inner boundary of a linear reservoir. It was not 
possible to lower the temperature in a step-wise 
manner with the method and equipment devised. 
Instead of a step-wise solution, therefore, a single 
pressure drop over one period of time will be used. 
Then, by introducing the reservoir diffusivity 
analog = k/C dt into the equation, it reduces 


to 
Q=rcgiP/ te... ... 2. 
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The dimensions of the model will be expanded to 
reservoir size by using a horizontal scale of 1 in.: 
1,000 ft and a vertical scale of 1 in.:100 ft. The 
cross-sectional area, numerically, then will be 
400,000 times the thickness of the model. The fluid 
was assigned a compressibility of 5.0 x 10-6 vol/ 
vol/psi and a viscosity of 0.4 cp. The rock and fluid 
properties then were related to the thermal properties 
by setting the conductivity, in Btu’s per (hour)(square 
feet)(degrees Fahrenheit per feet), equal to perme- 
ability in millidarcies and by setting the specific 
heat-density term equal to C p¢(10)®. The porosity 
is the only quantity which varies in the latter group 
of terms since each layer is assumed to contain the 
same fluid. 

For the purpose of this analogy, a reservoir time 
of 5,000 hours was chosen. Since no pressure change 
has been felt at the outer boundary of the reservoir 
during this time, the system may be considered in- 
finite. 

The copper plate is then analogous to a reservoir 
12,000-ft long, 4,000-ft wide and 25-ft thick, with a 
porosity of 25.5 per cent, a permeability of 229 md 
and the fluid properties previously mentioned. The 
initial reservoir pressure is considered to be well 
above the bubble point. 

The lead plate will represent a reservoir having 
the same horizontal dimensions and a pay thickness 
of 22 ft. Its porosity is 10.9 per cent and its per- 
meability 23 md. The fluid properties and initial 
reservoir pressure will be the same as those in the 
copper-plate analogy. 

To calculate the cumulative influx of fluid, a 
pressure drop of 100 psi was introduced suddenly 
at the interior boundary of each individual reservoir. 
The fluid produced from the reservoir simulated by 
copper would then be 








(2)(250(4000)(5 x 10-5)(255) —_ ‘iin 
5.61 7 


= 62,600 bbl. 


Similarly, the reservoir simulated by lead had a 
cumulative fluid influx of 11,320 bbl at the end of 
5,000 hours. Assuming that the two reservoirs were 
separated by a thin impermeable layer, the total 
fluid produced would be the sum of the two, or 73,920 
bbl. 

In calculating the fluid influx for this system 
without the impermeable layer between the reser- 
voirs, an average weighted porosity was used. The 
diffusivity was calculated by Eq. 1 and converted 
to a reservoir diffusivity analog. With all the 
other factors remaining the same, the cumulative 
fluid influx was 76,400 bbl. This is an increase of 
2,480 bbl, or 3.25 per cent, over the total cumula- 
tive production of the two layers calculated individ- 
ually. This increase could be attributed to the pres- 
ence of a cross flow between the two reservoirs. 

A third layer, represented by Wood’s metal, had 
thermal properties proportional to a permeability of 
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8.9 md, a porosity of 11.5 per cent and a diffusivity 
of 0.028 sq ft/sec. Its thickness is 12 ft and the 
other dimensions, as well as its fluid properties, are 
the same as for the copper and lead layers. Using 
the same pressure differential for the same 5,000 
hours, this reservoir would produce 3,930 bbl. When 
combined with the other two reservoirs to form a 
three-layered system susceptible to cross flow, the 
cumulative production from the system would be 
84,200 bbl. In comparison, the sum of the cumulative 
production from each layer would be 77,850 bbl. 
According to these calculations then, an increase 
in production of 6,350 bbl, or 8.2 per cent, is real- 
ized when the pays are not isolated by impermeable 
layers. 

Several other hypothetical stratified systems have 
been analyzed, and all have shown a higher cumula- 
tive fluid influx for the stratified systems where 
cross flow prevailed. 

Cross flow is apparent in the thermal model by 
the absence of a marked difference between the 
temperature distributions in the various plates of a 
given system. When plotted after the manner of 
Carslaw and Jaeger,’ the data points taken all fell 
on the same dimensionless time line, regardless of 
the material from which the data were taken. 

It seems reasonable then to consider the higher 
values, obtained using a common diffusivity, to be 
values which compensate for a cross flow. 

The effects of having a material of low conduc- 
tivity between plates of relatively high conductivity 
are illustrated in Fig. 8. The diffusivity measured 
on the aluminum side of the plate is high initially, 
falls rather rapidly and then begins leveling off as 
the steady-state condition is approached. The 
stainless-steel-side data indicate little change, if 
any, during the course of the run. 

The steep slope of the aluminum curve at low 
dimensionless times is due to a changing rate of 
cross flow induced by the large initial temperature 
differential between the rubber and aluminum. The 
steel, compared to aluminum, is itself fairly low in 
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FIG. 8 — EFFECT OF A LAYER OF LOW 
CONDUCTIVITY ON EFFECTIVE THERMAL 
DIFFUSIVITY. 


conductivity and experiences a rather small differ- 
ential across its thickness. 

From mathematical model studies of layered sys- 
tems, one finds that the rate of drainage may be 
expected to achieve a pseudo-steady state. The time 
required, rate of depletion and approach to steady 
state depends on the ratios of porosities, permea- 
bilities, layer thicknesses and lengths. 

To determine the effect of the thickness of the 
low-conductive layer on the performance of the sys- 
tem, runs were made on a brass, rubber and steel 
plate. The rubber thickness was increased from 
0.009 in. in the first set of runs to 0.106 in. in the 
second set of runs. These results are illustrated in 
Fig. 9(A and B). The data taken from the thin-layer 
system gave results similar to those obtained in 
two-layered systems but with a definite temperature 
differential across the rubber, i.e., a brass-rubber 
and a steel-rubber plate, but with the rubber acting 
as if it were higher in conductivity than it is. 
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FIG. 9— THE EFFECT OF VARYING THE THICKNESS 
OF A LOW-CONDUCTIVE MATERIAL BETWEEN TWO 
METALS ON DIFFUSIVITY. 
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The 0.106-in. rubber layer caused a higher initial 
diffusivity difference between the plates. These 
results are similar to those obtained with the 
aluminum-rubber-steel system. 


In Fig. 10(A and B) the temperature differential 
across the thickness of the plate is plotted, in di- 
mensionless form, with time. The maximum differ- 
ence is shown at early times close to the source. 
Peaks are indicated at X values of 0.8 and 0.6, 
while none are discernible at distances farther than 
this from the source. 

The effect of the thicker rubber layer is fairly 
obvious at early times. If the system were to be 
considered as a stratified pay with a maximum 
terminal-pressure differential of 200 psi, the peak 
at 0.165 dimensionless pressure would correspond 
to a pressure differential of 33 psi across the width 
of the sand. This differential would be felt at a 
distance of 200 ft from the pressure sink if the 
system were considered to be 1,000-ft long. At the 
same point and at the same time, the differential 
would be 10 psi in the thin-layer system. This dif- 
ference is less pronounced at later times. 

Using the same analogy at later times, the dif- 
ferential across the thick tight zone would be 7- to 
12-psi while across the thin tight zone the maximum 
differential would be 11l-psi near the source and 2- 
to 6-psi elsewhere. 


CONCLUSIONS 


Experimental confirmation has been obtained of 
the utility of an average permeability, porosity and 
compressibility for layered or stratified reservoirs 
in which cross flow occurs. This proper value is 
obtained by Eq. 5 

" 
2 knhn 


je et 
zs Ci dh by 
n=1 


A knowledge of the proper method of weighting 
core-analysis data has considerable application. In 
its simplest application in petroleum reservoir equa- 
tions involving transient pressure where rock com- 
pressibility is neglected and the rock is filled with 
a single fluid of constant compressibility, the proper 
single value of permeability to use in a stratified 
pay is indicated to be given by the sum of the per- 
meabilities at each 1-ft interval divided by the total 
pay thickness in feet. If data are available at ir- 
regular spacing, then the proper value of k is given 
by (Zkb) + Lh. The proper value of porosity to be 
used in the transient equation is equal to (2 ¢4)/ 
(Xb). 

Using these definitions, a layered reservoir with 
cross flow can be considered as a uniform homoge- 
neous reservoir of the same thickness and length if 
the permeability -porosity ratio is represented by 
(Lk b)/X db. 

It should be clearly emphasized that this recom- 
mendation applies when the reservoir approaches 
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pseudo-steady state because the effective value of 
k will be larger during the early transient period. 

In their simplest application, these findings in- 
dicate that a first estimate of the proper value of 
permeability and porosity may be obtained quite 
easily by preparing a table listing the permeability 
and porosity at each 1-ft interval. Sum each column 
and divide by the total pay thickness to obtain the 
average permeability and porosity. This practice 
has been almost universal and may represent a first 
approximation. 

In reality, an oil-water interface usually exists 
and the connate water varies from foot to foot. Since 
this is true, it should be recognized that the e/fec- 
tive compressibilities are different at each 1-ft in- 
terval. For this case, the proper porosity- 
compressibility product is obtained from (= c ¢ h)/ 
(Xb), where the compressibility c at each 1-ft in- 
terval represents the properly weighted value for 
the particular rock properties, porosity, connate 
water and oil. Although this calculation takes 
longer, it yields a better answer for the best single 
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FIG. 10 —— THE CHANGE IN TEMPERATURE DIFFER- 
ENTIAL ACROSS A LAYERED SYSTEM. 





numbers to be used to represent the effective prop- 


erties of the pay. 
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ABSTRACT 


Analytical solutions are obtained for calculating 
the pressure distribution in rectangular fields due 
to injection and/or producing wells located any- 
where within the field. The field is assumed to be 
homogeneous with either constant pressure or no- 
flow boundaries. These solutions are extended to 
include the calculation of average reservoir pres- 
sure and permeability from pressure build-up curves. 
Numerical examples are given to illustrate the 
application of equations to practical problems. 


INTRODUCTION 


Pressure distribution analyses are of consider- 
able importance in the field of reservoir mechanics. 
A number of papers dealing with this subject have 
appeared in the petroleum literature. Perrine! in 
1956 presented a summary of the methods available 
for pressure build-up calculations and discussed 
the applicability of each in some detail. More 
recently (1958), Hazebroek, et al,2 discussed a 
theoretical method for obtaining pressure distribu- 
tion in a radial field. Also in 1958, Nisle? gave a 
theoretical solution of the pressure distribution in 
a field extending to infinity with a partially pene- 
trating line source at the origin. 

Mathews, et al,4 applied the method of images to 
existing solutions for pressure distribution in radial 
fields and obtained expressions for calculating 
pressures in bounded reservoirs. However, a litera- 
ture survey revealed that a general analytical 
solution for the pressure distribution in a rectangular 
field, with injection and/or production wells located 
anywhere within the field, was not available. Such 
a solution could be used to approximate the pres- 
sure distribution in a field the shape of which is 
nearly rectangular. Furthermore, it could be applied 
to pressure calculations pertaining to a single well 
taken from a system of wells which were drilled in 
a rectangular pattern. In this case the boundaries 
of the well drainage area are considered to be the 
perpendicular bisectors of the lines joining the 
well and its nearest neighbors. 
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The basic problem considered in this paper con- 
sists of a rectangular field with either constant- 
pressure or no-flow boundaries and with either an 
injection or production well located anywhere within 
the field. The physical properties of the rock and 
fluids present are considered to be constant. Ana- 
lytical expressions are obtained for the pressure 
distribution within the field during the production 
or injection periods. These expressions, which are 
infinite series, can be used to evaluate pressure at 
any point within the field as a function of time and 
position. For fields containing more than one well 
(production and/or injection), the solutions for each 
well acting independently are superposed. That is, 
at any point in the field, the pressures due to each 
injection or production well acting alone can be 
algebraically added to give the pressure resulting 
from all wells acting simultaneously. Evaluation of 
the series solutions given in this paper is readily 
accomplished by means of an electronic digital 
computer; and, in most cases, sufficient accuracy 
is obtained by setting the upper limits of summation 
at 20. 

These solutions are extended to include the 
calculation of average reservoir pressure and per- 
meability from build-up data. A number of numerical 
examples are given to illustrate application of the 
solutions to actual field problems. 


ANALYTICAL SOLUTIONS 


The partial differential equation representing the 
pressure distribution in a homogeneous rectangular 
field (see Fig. 1), having a single source or sink 
of strength Q, can be written* 


2 2 
te =a SP + pQ(l.g).- ~~ (D 
ax2 dy2 t 


where 


B 
a = 157.952 2" and p = 686,905 —2* 


x, y = space co-ordinates, ft, 

(1, q) = location of source or sink, ft, 
p = pressure at (x, y), psi, 
@ = porosity, fractional, 





*This equation is the same as the two-dimensional flow 
equation with a source (see Ref. 5). 








compressibility of fluids plus rock struc- 
ture, psi~}, 


pt = viscosity, cp, 
k = permeability, md, 
t = time, days, 
Bo = oil formation volume factor, reservoir bbl/ 
STB, and 
QO = specific production (+) or injection (-) rate, 


STB/D-ft of thickness. 

Although Eq. 1 is for a homogeneous field with 
single-phase flow, by using ‘‘average’’ field prop- 
erties,© the solutions of this equation can be 
applied to multiphase-flow systems. 

Since Eq. 1 is a linear partial differential equa- 
tion, the principle of superposition of solutions 
applies. Consequently, if there are two or more 
production and/or injection wells, the pressures 
due to each well acting independently can be added 
algebraically to give the resultant pressure. 


CASE I — CONSTANT PRESSURE BOUNDARIES (CPB) 

The solution of Eq. 1, subject to the constant 
pressure boundary condition and an initial constant 
pressure distribution of P,, is given by the equation 
(see Appendix), 


p(x,y,t) = Po 
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where P, = uniform initial pressure distribution, 
psi, 
m,n = summation indices, integers, and 
a, b = dimensions of the rectangular section, 


ft. 


The infinite series of Eq. 2 represents the pres- 
sure distribution as a function of location (x, y) and 
time (t). Note that the exponential expression within 
the brackets of Eq. 2 is the time-dependent part of 
each term. For large values of time ¢, this expression 
will go to zero and the steady-state solution will 
result. The series converges rapidly, except at the 
location of the source or sink, where it diverges to 
infinity. 

Fig. 2 depicts the pressure distribution in a 
square field with an injection well at the center 
and an initial pressure distribution P, of zero. This 
figure is obtained by evaluating Eq. 2 along lines 
of constant ordinate. Note that the pressures in- 
crease as one moves toward the injection well and 
that the rate of pressure change increases similarly. 
Note also that the pressures at the boundaries re- 
main constant at zero, in accord with the boundary 
conditions. 
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CASE II — NO-FLOW BOUNDARIES (NFB) 

The solution of Eq. 1, subject to the no-flow 
boundary condition and an initial constant pressure 
distribution of P,, is as follows (see Appendix), 


p(zy,t) = Pp PQ [ts 2 F 
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. (3) 


Except at the location of the line source (where 
it diverges to infinity), Eq. 3 can be used to evalu- 
ate the pressure for any specified location at any 
given time. The series converges at a reasonable 
rate, and good approximations to the pressure can 
be obtained by setting the upper limits of the sum- 
mation to 20. 


Inspection of Eq. 3 shows that (in accord with 
the assumed initial condition) p(x, y, t) equals P, 
when ¢ is zero. It can also be seen that the term t/a 
increases with time while all other terms get smaller 
and smaller. Consequently, no steady-state condition 
is ever reached because the pressure at any location 
must always increase with time so long as constant 
injection rate is maintained. It is interesting to 
note that in Eqs. 2 and 3 the exponential expressions 
depend on time but not on location, whereas the 
trigonometric expressions depend on location but 
not on time. 


Fig. 3 depicts the pressure distribution in a square 
field with an injection well at the center and an 
initial pressure distribution P, of zero. Since there 
is no flow across the outer boundary, pressures in 
the field increase with time. After an initial transient 
period, a semisteady-state period is encountered 
whereby with increasing time the curves maintain 
the same shape but move upward along the pressure 
axis. Note that the slope of all curves at the bound- 
aries is zero, in accord with the no-flow boundary 
condition, 
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BUILD-UP ANALYSIS 


The solutions given in the previous section 
express the pressure distribution in the reservoir 
during the periods of production or injection. Once 
the production or injection terminates, they no longer 
express the pressure behavior of the system. To 
utilize these solutions after the termination of pro- 
duction or injection, we resort to the convolution 
integral,” and assume a period of constant produc- 
tion T prior to a shut-in period of t,. The following 
results. 


Ap = p(T + t,)—p(t,). dca esi 


where Ap = P, - ps (ts), 
bs = pressure after shut-in, psi, 


! 


p = flowing pressure as given by Eqs. 2 and 
3, psi, 

t, = build-up time, days, and 

T = production or injection time prior to 


shut-in, days. 
If the well drainage area is square (a = 5), the 
following dimensionless quantities® may be defined. 


= kAp 5) 
Ap = 0.001127 one 
QuB 


kt, 
0.006331 ——> 
pcpa 


Ca 
u 


. (6) 





k 
0.006331 — u 


T 
gcpua2 


Furthermore, if the pressure history of a well can 
be suitably expressed by either Eq. 2 or 3, then Eq. 
4 may be evaluated in terms of these dimensionless 
quantities. Fig. 4 summarizes the results of a number 
of such evaluations for a square section with either 
an off-center or a centrally located well. The bound- 
aries of the section are either at constant pressure 
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FIG. 1 — RECTANGULAR FIELD. 
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or have no flow. In these calculations a point 
corresponding to a wellbore radius of 3-in. was used. 
It was observed, however, that varying this radius 
from 3 to 10 in. had no significant effect on the 
evaluation of Ap. 

Fig. 4 includes the dimensionless pressure-vs- 
time plots for a production well located at / = q =4. 
It is noted that the slope of these curves and the 
curves of centrally located wells are the same. It 
is further observed that the curve of NFB case is 
higher than that of a centrally located well. This 
will result in lower pressures during the build-up 
period. This is reasonable because during the pro- 
duction period the effect of boundaries is greater 
at the well (compared to centrally located wells), 
resulting in lower flowing pressures at the time of 
shut-in. This, in turn, results in lower build-up 
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FIG. 2—PRESSURE DISTRIBUTION FOR FIELDS WITH 
CONSTANT PRESSURE BOUNDARIES. 
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FIG. 3— PRESSURE DISTRIBUTION FOR FIELDS WITH 
NO-FLOW BOUNDARIES. 


225 




















CONSTANT PRESSURE BOUNDARIES (CPB) 
—--- NO-FLOW BOUNDARIES (NFB) 


T=DURATION OF PRODUCTION 











KAP 
QuBo 





AP=0.001127 





kt 


T= 0.00633! Sanat 


FIG. 4 — DIMENSIONLESS PRESSURE BUILD-UP 
CURVES. 
pressures. The curve of CPB case of Fig. 4 is 
lower than that of a centrally located production 
well. This will result in higher build-up pressures. 
This is also reasonable since during the production 
period the constant pressure boundaries (at the 
original field pressure) are closer to the producing 
well, resulting in higher flowing pressures at the 
time of shut-in. This higher flowing pressure will 
result in higher pressures during the build-up period. 
For t < 0.02, observe that all curves of Fig. 4 
have the same constant slope; this slope is 


A(Ap) 

A(logt) 
Since in practice the actual pressure against time 
is plotted on a semilogarithmic scale, we equate 


this slope (with proper conversion factors) to that 
of Eq. 7. 


TENE x se ee 1 


ho, 
QuB, |A(logts)| FIELD 


. (8) 





0.1675=0.001127 


From Eq. 8 the effective permeability k is 
k = 148.62 vie, 
; [p/A(logt,) rreun’ * 


With the aid of Fig. 4 and Eqs. 5, 6 and 9,* a 
number of calculations can be made utilizing the 
pressure build-up curves of practical field cases. 
A summary of possible calculations is given in 
Table 1. To illustrate the application of the method, 
a numerical example is worked out for each case. 
In all of the following examples the field is assumed 








*Comparing Eq. 9 above and Eq. 6 of Ref. § one notes a 
discrepancy of about 10 per cent, This is due to the effect of 
boundaries and also the fact that in the curves of the reference 
a steady-state condition is assumed prior to shut-in. This is 


equivalent to setting P(T + t,) = P(T) in Eq. 4. 


226 





TABLE 1 — SUMMARY OF POSSIBLE CALCULATIONS 
FOR FOUR CASES 











Case | Case 2 Case 3 Case 4 
GIVEN: P,,k, d, c Buildeup Build-up Complete 
by Q, Bo, a Curves Curves ae 
LQ, By d, Git, + yee 
Q, Bo. = rag H, 
1 Bo 
FIND: Build-up k k, Field k, a 
Curves Pressure 





to be square with a production well at the center. 


CASE 1 

Given: PR, = 3,500 psi, k = 1.5 md, 6 = 0.20, c = 
1.5 x 10~5psi-}, p = 0.33 cp, Q = 1.0 STB/ 
D-ft, Bg = 1.5 reservoir bbI/STB and a = 
2,000 ft. Assume steady state before shut-in 
for CPB and a production time of 365 days 
for NFB. 

Find: the pressure build-up curves. 

The production time of NFB, from Eq. 6, corres- 
ponds to T = 0.88. Suppose we want to obtain build- 
up pressures at ee se days. From 
Eq. 6, ts = 10-2 days corresponds to 


1.5 x 1072 
0.20 x 1.5 x 1075 x 0.33 x (2,000)2 


= 2.4 x 107°. 
From Fig. 4 this value corresponds to 
Ap = 0.63 CPB 
and _ 
Ap = 1.42 NFB. 
Using Eq. 5, Ap corresponding to the above will be 
Ap = 185 psi CPB 
and 
Ap = 415 psi NFB, 
or the actual pressures are 
Ps = 3,500 - 185 = 3,315 psi CPB, 
and 
bs = 3,500 — 415 = 3,085 psi NFB. 
The values of pressures can be obtained in a similar 
manner for the rest of the ¢, values. A plot of the 
build-up curve is shown on Fig. 5. 





t = 0.006331 


CASE 2 
Given: p = 0.33 cp, Q = 1.0 STB/D-ft, B, = 1.5 
reservoir bbl/STB and the pressure build-up 
curves of Fig. 5. 
Find: the effective permeability &. 
The slope of the build-up curves of Fig. 5, taking 
the section between t =10 and 1, is 


[ Ap ] - 3,408 ~ 3,360 _ 4p. 
A (log &) Jeretp O~(~2) 


From Eq. 9 we have 





1.0 x 0.33 x 1.5 
48 





k = 148.62 = 1.53 md. 


CASE 3 
Given: @ = 0.20, c = 1.5 x 1075 psi~!, a = 2,000 ft 
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FIG. 5— PRESSURE BUILD-UP CURVES OF THE 
FIELD PROBLEM. 
and the data of Case 2. 
Find: the effective permeability & and field pres- 
sure p. 

The value of k is obtained by the procedure given 
in Case 2. Let us obtain the values of pressure at 
time t = 10~! days from Fig. 5. 

pb; = 3,360 psi CPB, 
and 

b2 = 3,137 psi NFB. 

The dimensionless time corresponding to this ¢ = 
10-1, using the value of k calculated in Case 2, is 


1.53 x 107! 
0.20 x 1.5 x 107° x 0.33 x (2,000)? 


= 2.42 x 1074. 


From Fig. 4 this corresponds to Ap of 0.47 for CPB 
and Ap of 0.37 for NFB. Note that these Ap values 
are the differences of Ap values corresponding to 
t = 2.42 x 1074 and the Ap value of steady- or 
semisteady-state condition on the same curve. From 
Fig. 4 the steady-state Ap for CPB is zero, while 
for NFB it ranges between 0.10 and 1.0 depending 
on the duration of production. Using these Ap val- 
ues from Eq. 5, the average stabilized field pres- 
sures are 


t = 0.006331 








p = 3,360 + 0.47 x 1.0 x 0.33 x 1.5 
1.53 x 0.001127 
= 3,496 psi CPB, 
and 
p = 3,137 + 9:32 x1.0 x 0.55 x15 
1.53 x 0.001127 
= 3,244 psi NFB. 


CASE 4 

Given: @ = 0.20, c = 1.5 x 1075 psi=!, p = 0.33 ep, 
Q =1.0 STB/D-ft, B, = 1.5 reservoir bbl/ 
STB and the complete pressure build-up 
curves of Fig. 5. 

Find: the average permeability k and the dimension 
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a of the square field. 

The average permeability k is obtained as in 
Case 2. From Fig. 5, the pressure increase from 
t, = 10-1 (transient period) to t, = 4 x 102 (steady- 
State period) is found to be 

Ap = 3,496 — 3,358 =138 psi or Ap =0.485 for CPB, 
and 

Ap = 3,244 —3,126 =118 psi or Ap =0.360 for NFB. 
In Fig. 4, these values of Ap are seen to correspond 
to 

F = 2.4 x 10~* for CPB, 
and 

F = 2.0 x 10~ for NFB. 
Solving Eq. 6 for a and using t, = 10~!, we obtain 
a = 2,200 ft for NFB and a = 2,000 ft for CPB. The 
difference in the a values computed here is due 
principally to errors introduced when reading the 
graphs. 


SUMMARY 


Analytical solutions are obtained for calculating 
the pressure distribution in rectangular fields due 
to injection and/or production wells located any- 
where within the fields. Boundaries of the field are 
assumed to be either at constant pressure or have 
no-flow. Properties of the field and of the fluids 
present are considered to be constant. The solu- 
tions are extended to include the calculation of 
average reservoir pressure and permeability from 
pressure build-up curves. 

It is observed that, in most of the literature on 
pressure distribution analysis, the well is located 
at the center of the field. The solutions given in 
this paper are perfectly general and apply to fields 
with either off-center or centrally located wells. 
The analysis presented in this paper can also be 
used for well-interference problems. 


NOMENCLATURE 


a, b = sides of the rectangular section, ft 
B, = oil formation volume factor, reservoir bbl/ 
STB 
c = compressibility of rock and fluids present, 
psi 
CPB = constant pressure boundaries 
k = effective field permeability, md 
l, q = location of source or sink Q, ft 
m,n = summation indices, integers 
NFB = no-flow boundaries 
P, = initial pressure throughout the field, psi 
ps = pressure after shut-in, psi 
p = pressure at location (x, y), psi 
Ap = pressure during build-up, psi 
Ap = dimensionless pressure, Eq. 5 
Q = specific production (+) or injection (-—), STB 
/D/ft of thickness 
= time, days 
dimensionless time, Eq. 6 
build-up time, days 


t 

t 
ts 
T 


if 


duration of production or injection prior to 
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shut-in, days 

dimensionless, duration of production or in- 
jection prior to shut-in, Eq. 6 

space co-ordinates, ft 


157.952 2c 


886.905 Hot 


= porosity, fraction 
viscosity, cp 
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APPENDIX 


We define the symbols T’ and Q’ as follows. 
t 


r ag, =. 


Using these symbols in Eq. 1 and representing the 
point source Q by the Dirac delta function (5),” we 
have 


2 2 > 

Lal r <i = OP + Q'8(x-1)8(y-9). (10) 
2 dy2 oT 

The boundary and initial conditions of Eq. 10 are 
p(0,y,T’) p(a,y,T’) = 0, 


p(x,0,T’) = p(x,b,T’) = 0, 


p(x,y,0) 


OK<y <b; 
0 <x «a; 
= Po(x,y), for all x and y. 
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Multiplying both sides of Eq. 10 by 


‘ mux , nwy 
sin“ sin ir 


a 
and integrating between the limits indicated, we 
have 


ff 


a 
f =, sin sin 
0 3 a 


re) 


d2p mux nw 
+ ~t) sin — sin i dxdy = 
dy2 a 


p 
dx2 b 


 dedy + Q' 
nwy 


sin big dxdy 
(11) 


5(x-l)8(y-q) sin 
( a 
Applying the finite Fourier sine transform 9 


o ; ba 4 «Max 
P,(m,n,T’) = [f, P(zey.T ; 


. nwy 
sin cs Gear. ou « 


to Eq. 11 and utilizing the given boundary condi- 
tions, we have 


(== 5 ( Tr d _ 
-7T —_ _—_— = ene 
at 6-62 Rater 2 dT’ Ps 


(m,n,T’) + Q’sin = sin =» A 
The solution of Eq. 13 lines to the —— initial 
condition can be written 
Q’sin mal/a sin naq/b 

w2 (m2 /a2 + n? /b2 ) 





P,(m,n,T)= es 


ba 
. max  — nwy 
+ff Po(x,y)sin ; a. ards | 


2 2 


i=. «oie 
x exp | -T a aay 





4 Q’sin mul/a sin ntqg/b =) 
st(at/at + ni/b*) 

The inverse of this transformed equation can be 
obtained from the Fourier Integral theorem. This 
transformation, subject to a constant initial pressure 
distribution P,, results in Eq, 2. 

The NFB solution, given by Eq. 3, is obtained 
by using the following boundary and initial condi- 
tions: 


2 5(0,y,T')= > p(a,y,T’)=0, 08 y <b5 
ox dx 


© 5(x,0,T’)=2- p(x,b,T’ )=0, Oex <a; 
dy dy 


p(x,y,0)=Po(x,y), for all x and y 


and the Fourier cosine transform. This procedure 
is similar to the one given above. aad 


SOCIETY OF PETROLEUM ENGINEERS JOURNAL 
























Deflocculation of Fractionated Montmorillonite 


By Sodium Polyphosphates 


F. W. JESSEN 
MEMBER AIME 
FARUK H. TURAN 


ABSTRACT 


The gel strength and viscosity of two different 
suspensions of fractionated montmorillonite clay 
were measured by using a Stormer viscosimeter 
and Fann V-G meter. The amount of sodium poly- 
phosphate (sodium hexametaphosphate, sodium 
tetraphosphate and sodium acid pyrophosphate) 
adsorbed by the clay particles was determined. 

A relationship seems to exist between the amount 
of polyphosphate adsorbed and the nature of the 
surface, edge or plate, assumed to be the principal 
sites for selective adsorption. Reduction in vis- 
cosity and gel strength was found to be directly 
related to the adsorption. Both chemical bonding 
and physical adsorption seem involved in defloccula- 
tion. Some indication of the axial ratios of particles 
of .1l- to .04-micron average size is possible 
through study of the adsorption mechanism. 


INTRODUCTION AND PURPOSE 


A knowledge of surface areas of clay particles 
in a suspension and a precise determination of 
phosphate adsorption should provide a basis for 
a better understanding of reduction in viscosity 
and gel strength caused by molecularly dehydrated 
phosphates. 

The purpose of this investigation was to study 
the adsorption of sodium polyphosphates on sus- 
pensions of fractionated sodium montmorillonite. 

Although there is considerable uncertainty re- 
garding details of the structure of montmorillonite, 
a generally accepted one is that suggested by 
Hofmann, Endell and Wilm.! The structure shown 
in Fig. 1 is an idealized one, and Marshall? has 
postulated that it is possible for several isomorphous 
replacements to occur. Thus, aluminum may replace 
silicon, in which case an exchangeable cation 


, would go along with the aluminum to neutralize 
what would otherwise be a net negative charge on 


the structure. 





Original manuscript received in Society of Petroleum 
Engineers office May 27, 1960. Revised manuscript received 
Sept. 18, 1961. 
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THE U. OF TEXAS 
AUSTIN, TEX. 


There is disagreement among various investigators 
as to the particle size and shape of montmorillonite 
clays. Marshall? reported the following results for 
the mean lateral dimensions of a Putnam clay 
(montmorillonite) — 3,730 A°” for particles having 
average diameters ranging from 500 to 1,000 A°, 
and 750 A° for those with an average diameter less 
than 200 A®°. The latter particles were estimated at 
50A° in thickness, which would correspond to a 
stack of five unit cell layers. On the other hand, 
Shaw3 has interpreted micrographs as indicating 
that particles of montmorillonite are readily broken 
down to the individual cell layers by dispersion in 
water. Brindley4 and Mering® indicate the ultimate 
particles to be of the order of 300 A®° in mean 
lateral dimension. Van Olphen® suggested the thick- 
ness of particles at four to five unit cell layers. It 
was suggested also that this number tends to become 
smaller with increasing dilution. Recently, Kahn 
and Lewis? determined the size of sodium mont- 
morillonite particles in aqueous suspension from 
electro-optical bi-refringence decay curves. The 





* Angstrom unit. 
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FIG. 1 — MONTMORILLONITE STRUCTURE SUG- 
GESTED BY HOFMANN, ENDELL AND WILM (REF. 1). 


229 


results were interpreted in terms of the diameter of 
an oblate spheroid, the lateral dimension of the 
smallest particle being 5,000 A°. Melrose 8 concluded 
that, although numerous electron micrographs of 
montmorillonite clay have been published, in none 
of these can any particles of definite outline be 
identified or measured. 


LABORATORY INVESTIGATION 


Viscosity measurements were obtained with the 
Stormer and Fann V-G instruments, the latter also 
serving to give the gel strength. A standard API 
filter press was used to obtain samples of filtrate 
from the suspension, which were further clarified 
by filtering through Whatman analytical-grade filter 
paper. 

A Lumerton photoelectric colorimeter (Model 400- 
A) was used in the colorimetric determination of 
phosphate. Known concentrations of sodium acid 
phosphate (NaH, PO,) representing 0.125, 0.250, 
0.625, 1.250, 2.500 and 5.000 ppm of phosphorus 
were used for calibration. In the determination of 
phosphate in each filtrate sample, three different 
readings were made. In no instance was the variation 
in individual determinations of the percentage 
transmission greater than 1 per cent in the range 
of 90 per cent transmission, while a 4 per cent vari- 
ation was found in the range of 50 per cent trans- 





TABLE 1 — 6 PER CENT SUSPENSION OF M FRACTION 
TREATED WITH SODIUM HEXAMETAPHOSPHATE 


Stormer Plastic 
Viscosity Viscosity 


(Ib/bb!) (Ib/bbl) (cp) (cp) 


Amount Amount 
Added Adsorbed Yield Point 


(Ib/ 100 f+?) 





0. - 46 32 
0.10 0.075 20 26 
0.25 0.182 15 22 
0.50 0.325 13 20 
1.00 0.556 10 20 
2.00 0.890 10 20 
4.00 1.120 14 22 


mission. Dilution of the filtrate was controlled so 
that determinations were made with 85 to 95 per cent 
transmittance, assuring an accuracy of 1 per cent 
in the phosphorus determination. 

Two clay suspensions were prepared using frac- 
tionated portions of a Wyoming bentonite studied 
by Jessen and Mungan.? The average particle diam- 
eter of one fraction (designated M) was 0.04 micron, 
while the other fraction (designated L) has an 
average particle size of 0.10 micron. 

Quantities equivalent to 0.10, 0.25, 0.50, 1.00, 
2.00 and 4.00 lb/bbl (1 1b/bbl = 2.8525 gm/liter) of 
sodium hexametaphosphate were added to the clay 
suspensions. The suspensions were stirred thor- 
oughly for 30 minutes, and the Stormer viscosity, 
the plastic viscosity and the yield point were 
determined. 

The filtrate was diluted with water and sulfuric 
acid added. The sample was allowed to boil in an 
Erlenmeyer flask for at least two hours to convert 
the metaphosphate to orthophosphate. !° An aliquot 
was used in the colorimetric determination of phos- 
phate, and the amount of the meta-, tetra-, and 
pyrophosphate adsorbed by the clay particles was 
calculated by subtracting the amount in the filtrate 
from the original quantity present. 

The results are presented in Tables 1 through 6 
to show the amount of sodium hexametaphosphate, 
sodium tetraphosphate and sodium acid pyrophos- 
phate adsorbed vs the change in yield point and 
the reduction in plastic viscosity and Stormer 
viscosity. 

Tables 7 through 12 present the extent of adsorp- 
tion of sodium hexametaphosphate on the clay 
particles calculated on the basis of varying axial 
ratios (cba). 


RESULTS AND DISCUSSION 


The interaction of polyphosphates with mont- 
morillonite has been studied by Van Olphen,!! whose 





TABLE 2 — 7 PER CENT SUSPENSION OF L FRACTION 
TREATED WITH SODIUM HEXAME TAPHOSPHATE 


Stormer Plastic 
Viscosity Viscosity Yield Point 
(Ib/bb!) (Ib/bbi) (cp) (cp) (Ib/ 100 +2) 
0. - 57.0 38 65 
0.10 0.078 34.0 28 29 
0.25 0.196 23.0 26 19 
0.50 0.320 18.0 25 12 
1.00 0.566 17.5 24 10 
2.00 0.900 14.5 22 9 
4.00 1.160 14.0 21 7 


Amount Amount 
Added Adsorbed 





TABLE 4 — 7 PER CENT SUSPENSION OF L FRACTION 
TREATED WITH SODIUM TETRAPHOSPHATE 


Plastic 
Viscosity 


Stormer 
Viscosity 
(Ib/bb1) (cp) 


Amount Amount 
Added Adsorbed 
(Ib/bbl) 
0 =~ 
0.10 0.078 
0.25 0.182 
0.50 0.345 
1.00 0.612 
2.00 1.110 
4.00 1.780 


Yield Point 
(1b/100 f+?) 











TABLE 3 — 6 PER CENT SUSPENSION OF M FRACTION 
TREATED WITH SODIUM TETRAPHOSPHATE 


Amount Amount Stormer Plastic 
Added Adsorbed Viscosity Viscosity 
(Ib/bbl) (Ib/bb1) (cp) (cp) 


Yield Point 
(Ib/100 ft2) 


TABLE 5 — 6 PER CENT SUSPENSION OF M FRACTION 
TREATED WITH SODIUM ACID PYROPHOSPHATE 


Stormer Plastic 
Viscosity Viscosity 
(Ib/bb!) (Ib/bbl) (cp) (cp) 


Amount Amount 
Added Adsorbed Yield Point 


(Ib/ 100 ft?) 





0. - 46 33 51 
0.10 0.078 26 Kt] 20 
0.25 0.179 18 24 17 
0.50 0.331 14 22 iB 
1.00 0.617 1 20 6 
2.00 1.000 1 20 6 
4.00 1.292 14 22 it 
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0 - 46 31 50 
0.10 0.082 20 26 9 
0.25 0.195 18 23 13 
0.50 0.350 14 22 10 
1.00 0.690 13 22 9 
2.00 1.176 13 22 8 
4.00 1.776 18 25 VW 
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findings favored physicochemical adsorption as an 
explanation for deflocculation. 

Sodium polyphosphates are linear polymeric 
structures composed of linked phosphorus-oxygen 
tetrahedra, with the general formula Na,,2P,03n41- 
page radii for tetrahedral coordination are 

.35 and O = 1.49 A°.!? Using a P-O length of 
i 84 A°, the base of a PO, tetrahedron is calculated 
to be 15 sq A°.'5 To compute the coverage of any 
specific clay surface with phosphate molecules, 
the axial distance (c, 6 and a) must be known. 
Melrose8 has suggested lateral dimensions (6 and 
a) in the ratio of 1:2 to 3:4 for clay particles. The 
thickness c of the smallest clay particles is not 
known exactly, but minimum values from 10 to 
50 A° have been suggested in the literature. These 
values define a particle of one to five unit cells 
in thickness. Assuming that fractionM with average 
particle diameter of 0.04 micron (400 A®°) is composed 
of particles five-unit-cells thick, the axial lengths 
are arrived at as follows. 

Volume of Particle = V = (7/6)d° = (3.14/6)(400)9 


= 33.75 x 10% cu A® 


Lateral Area =(b xa) =V/c = (33.75 x 106)/50 
= 67.5 x 104 sq A®°. 
Therefore, 
for b:a = 1:2, b = 580 A° and a = 1,160 A®, 
giving cb:a = 1:10:20; 
and 
for b:a = 3:4, b = 710 A° and a = 950 A®, 
giving cb:a = 1:14:19, 

In the first case, a mean lateral dimension of 
(580 + 1,160)/2 = 870 A® is indicated; and in the 
second case, (710 + 950)/2 = 830 A®. The axial 
ratios obtained here may be dssumed to apply to 





TABLE 6 — 7 PER CENT SUSPENSION OF L FRACTION 
TREATED WITH SODIUM ACID PYROPHOSPHATE 


Amount Amount Stormer Plastic 








Added Adsorbed Viscosity Viscosity Yield Point 

(Ib/bb!) (Ib/bb!) (cp) (cp) (1b/100 ft2) 
0 - 57 38 65 
0.10 .083 35 Kt) 30 
0.25 -194 22 2B 19 
0.50 -333 18 25 1 
1.00 667 13 22 8 
2.00 1. 166 10 20 6 
4.00 1.776 12 21 7 


particles of fractions with greater average sizes 
as well. Should particles be assumed to be of : 
single unit cell layer as postulated by Shaw,° 
resulting axial ratios (c:b:a) of 1:130:170 or 
1:100:200 will be obtained, and the mean lateral 
dimension is 1,650 A°. Similarly for the L fraction 
which has an average particle size of 0.10 micron, 
computations show 140 x 1,400 x 2,800 A° for the 
1:10:20 ratios, 120 x 1,680 x 2,280 A®° for 1:14:19 
ratios, and 35 x 3,500 x 7,000 A° for the 1:100:200 
ratios of c:b:a. Mean lateral dimensions range from 
2,000 to 5,200 A°, while the indicated thickness 
varies from 3 to 14 unit cell layers. All of these 
derived particle dimensions fall within the range 
of sizes reported by most investigators. It is apparent 
that a rather wide vafiation in adsorption might 
exist if the assumption is made that the principal 
action of polyphosphate molecules is on the broken 
edges. 

A compilation of the data is presented in Tables 
7 through 12. The coverage on any clay surface was 
calculated from the amount of clay fraction treated, 
corresponding to a definite number of particles with 
the afore-mentioned postulated dimensions, and the 
surface area of the polyphosphate adsorbed. Spe- 
cifically, consider 6 per cent M fraction treated 
with 0.10-lb/bbl polyphosphate (Table 2). 

Using as a basis 100 cc of suspension 


ec xb xa = 50 x 500 x 1,000 A®, 


Vol. of Clay 
Vol. of Particles 





No. Clay Particles = 


(60/2.3) x 1024 
33.75 x 106 





= 7.75 x 10!7, 


No. Polyphosphate Tetrahedra = 
(0.075) (2.8525) 
102 





x 6.02 x 1023 = 12.5 x 1020, 


where 2.8525 is conversion factor from pounds per 
barrel to grams per liter. 


Total Area of Polyphosphate = 12.5 x 1029 x 15 
= 18.8 x 107! sq A° 





TABLE 7 — NUMBER OF MOLECULAR LAYERS OF SODIUM HEXAMETAPHOSPHATE ON THE SURFACES 
OF 7 PER CENT L FRACTION WITH THE AXIAL RATIOS OF 1:14:19 


Number of Layers on the Surfaces of 








100 

Amnent Adbéed Amount Adsorbed Faces 
(Ib/bbl) Ib/bbI % A, (be) 
0. 10 0.078 78.0 0.760 
0.25 0.194 77.7 1.900 
0.50 0.320 64.0 3.130 
1.00 0.570 57.0 5.600 
2.00 0.900 45.0 8.820 
4,00 1.160 21.0 11,450 
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010 001 
Faces Faces 
Az (ac) Ag (ab) A, +A; A, +Az+Azg 
0.575 0.036 0.328 0.032 
1.435 0.090 0.818 0.081 
2.370 0.148 1.350 0.133 
4.220 0.264 2.410 0.238 
6.670 0.417 3.800 0.375 
8.600 0.537 4.900 4.830 

































































































































































































































Surface Area (ab) = 2 xa x b = 2 (500) (1000) 
= 106 sq A®. 

Total Surface Area (ab) = (7.75 x 1017) (10°) 

= 7.75 x 1023 sq A°. 


18.8 x 1021 _ 


Coverage = = 0.02. 
“ 7.75 x 1023 


For clay particles of 100-millimicron average 
size (L fraction), reduction in viscosity and yield 
point was obtained with all the polyphosphates 
employed at about equal concentrations. From 
Tables 7, 8 and 9, single-layer coverage by sodium 
hexametaphosphate adsorption on the edges of the 
L-fraction particles is found at 0.24, 0.27, and 0.13 
lb/bbl for corresponding particle dimensions having 
axial ratios of 1:14:19, 1:10:20 and 1:100:200, 
respectively. For the M-fraction particles, 0.52, 
0.57 and 0.27 lb/bbl of sodium hexametaphosphate 
are required to cover the edges of particles with 


like axial ratios (Tables 10, 11 and 12). Substantially 
complete deflocculation of the clay is observed 
with coverage of the edges, as evidenced by the 
destruction of gel properties and the low viscosity 
values. 

The experiments of Michaels!3 showed that the 
most likely site for adsorption of polyphosphates 
by kaolinite is the exposed edge of the gibbsite 
layer of a fractured sheet. Should the adsorption 
of polyphosphates on montmorillonite occur in the 
Same manner, some degree of selectivity may be 
expected, particularly if direct chemical bonding 
takes place. It follows, therefore, that sharp 
changes in viscosity and gel properties would 
result because the neutralization of positively 
charged sites by attraction of anions would sub- 
stantially reduce edge-to-face, as well as edge- 
to-edge, configurational positions. 

Van Olphen!! found the order of magnitude of 
adsorption of sodium hexametaphosphate to be about 
14 m. eq.*/ 100 gm of unfractionated sodium montmoril- 
lonite. This corresponds to 0.35 lb/bbl of 7 per 
cent clay suspension. 





* Abbreviation m. eq. is milliequivalent. 









TABLE 8 — NUMBER OF MOLECULAR LAYERS OF SODIUM HEXAMETAPHOSPHATE ON THE SURFACES 
OF 7 PER CENT L FRACTION WITH THE AXIAL RATIOS OF 1:10:20 


Number of Layers on the Surfaces of 











Aument dbded Amount Adsorbed sal 
(Ib/bbl) Ib/bb! % A, (be) 
0.10 0.078 78.0 0.884 
0.25 0.194 77.7 2.210 
0.50 0.320 64.0 3.650 

1.00 0.570 57.0 6.500 
2.00 0.900 45.0 10.200 
4.00 1.160 21.0 13.200 








010 001 

Faces Faces 

Az (ac) Ag (ab) A, +A, A, +Az+Az3 
0.442 0.0442 0.295 0.038 
1.105 0.0111 0.737 0.097 
1.825 0.1835 1.215 0.160 
3.250 0.3250 2. 160 0.284 
5.100 0.5100 3.400 0.448 


6.600 0.6600 4,400 





0.575 











TABLE 9 — NUMBER OF MOLECULAR LAYERS OF SODIUM HEXAMETAPHOSPHATE ON THE SURFACES 
OF 7 PER CENT L FRACTION WITH THE AXIAL RATIOS OF 1: 100: 200 


Number of Layers on the Surfaces of 











100 

Amount Added Amount Adsorbed Faces 
(Ib/bbl) Ib/bbI % A, (be) 
0.10 0.078 78.0 1.900 
0.25 0.194 77.7 4.750 
0.50 0.320 64.0 7.850 
1.00 0.570 57.0 10.400 

2.00 0.900 45.0 22.000 
4.00 1. 160 21.0 28.500 








010 001 
Faces Faces 
Az (ac) Ag (ab) A, +Az A, +Az+Azs 
0.950 0.010 0.633 0.001 
2.375 0.024 1.580 0.024 
3.925 0.040 2.610 0.040 
5.200 0.052 4.660 0.051 
11.000 0.111 7.330 0.110 
14.250 0.143 9.470 0.142 





TABLE 10 — NUMBER OF MOLECULAR LAYERS OF SODIUM HEXAMETAPHOSPHATE ON THE SURFACES 
OF 6 PER CENT M FRACTION WITH THE AXIAL RATIOS OF 1:14:19 







Number of Layers on the Surfaces of 








dumeent Cdided Amount Adsorbed ral 
(ib/bbl) Ib/bbl Fe A, (be) 
0.10 075 75.0 0.345 
0.25 0.183 73.0 0.812 
0.50 0.325 65.0 1.440 
1.00 0.556 55.6 2.470 
2.00 0.890 44.5 3.960 
4.00 1.120 28.0 4.980 






010 001 

Faces Faces 
Az (ac) Ag (ab) A, +A2 A, +Az+Az3 
0.258 0.018 0.148 0.016 
0.608 0.043 0.348 0.038 
1.080 1.077 0.620 0.068 
1.850 0.132 1.060 0.117 
2.970 0.210 1.700 0.187 
3.730 0.264 2.140 0.235 
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The amounts of sodium hexametaphosphate re- 
quired to cover the broken edges for the 6 per cent 
M-fraction and 7 per cent L-fraction clays are 0.27 
and 0.57 lb/bbl, respectively, when the assumed 
particle axial ratio is 1:10:20. The average of the 
two values is 0.42 lb/bbl. If the axial ratio of 
clay particles is taken to be 1:14:19, an average 
of 0.38 lb/bbl of sodium hexametaphosphate should 
be adsorbed. This latter value is in good agreement 
with the value found by Van Olphen and gives 
support to the postulated specificity of adsorption 
on the edge surfaces. With axial ratios of 1:100:200, 
and correspondingly less broken edge surface, 
average adsorption for the two fractions is computed 
to be 0.20 Ib/bbl. 

Now, if the action be chemisorption, the phosphate 
tetrahedra would be attached to the aluminum by 
replacing the hydroxyl ions on the edges. There 
are two hydroxyl ions for each aluminum atom, the 
distance between the latter being 2.6 A° (Fig. 1). 
There are two aluminum atoms and four hydroxyl 
ions for every 7.8 A° unit repeating section; this, 
multiplied by the unit thickness of 14 A®°, is 109.2 
sq A°. Hence, the area can be covered completely 
by 109.2/15 = 7.3 phosphate tetrahedra. On the 
other hand, since there are only four replaceable 
hydroxyl ions in this same area, an area of only 
60 sq A° may be assumed to be covered; or, in other 
words, only 55 per cent of the area may thus be 
occupied. In the latter case, the amount of sodium 
hexametaphosphate required to replace all hydroxy] 
ions on the edges is 0.14, 0.15 and 0.07 lb/bbl for 
the 7 per cent L fraction, and 0.28, 0.31 and 0.15 
lb/bbl for the 6 per cent M fraction (as before, con- 
sidering axial ratios of 1:14:19, 1:10:20 and 1:100: 
200, respectively). 


If the principal selectivity of polyphosphate ions 
is by chemical bonding, then approximately 60 per 
cent reduction in the yield point of the L-fraction 
material is attained when this phase is reached; 
in the case of the M-fraction clay, greater edge 
surface area allows for more such direct bonding, 
and 80 per cent of the reduction in yield point may 
be attributed to such action. The data show that 
deflocculation results with complete coverage of 
the broken edges, and no further reduction in yield 
point is obtained regardless of the amount of poly- 
phosphate further added. For axial ratios of 1:100: 
200, saturation of the edge surface is complete when 
the corresponding reduction in yield point is 60 
per cent. Since substantial quantities of poly- 
phosphates still are adsorbed, this may be inter- 
preted to mean one of the following: (1) adsorption 
takes place on the plate (flat) surfaces — even 
though in these experiments with an excess of 
chemical, coverage of the .001 surface (ab) is com- 
puted to be only .143 molecular layers; (2) multi- 
molecular layers are adsorbed on the edge surfaces; 
or (3) the configuration for the 1:100:200 axial 
ratio is rather improbable. 

Sharper decreases in viscosity and gel strength 
are noted with sodium acid pyrophosphate. This 
may be due to comparatively better coverage of the 
edge sites since less interference may be occasioned 
in the attachment of the relatively smaller pyro- 
phosphate. 


SUMMARY 


The adsorption of sodium hexametaphosphate, 
sodium acid pyrophosphate and sodium tetraphos- 
phate on fractionated clay suspensions — primarily 
sodium montmorillonite — has been determined. 





TABLE 11 — NUMBER OF MOLECULAR LAYERS OF SODIUM HEXAMETAPHOSPHATE ON THE SURFACES 
OF 6 PER CENT M FRACTION WITH THE AXIAL RATIOS OF 1:10:20 


Number of Layers on the Surfaces of 





100 
Amount Adsorbed 





Amount Added 


Faces 
(Ib/bb!) Ib/bbI ba A, (be) 


0.10 075 75.0 0.408 
0.25 0.183 73.0 0.960 
0.50 0.325 65.0 1.710 
1.00 0.556 55.6 2.920 
2.00 0.890 44.5 4.680 
4.00 1,120 28.0 5.900 


010 001 
Faces Faces 
Az (ac) A; (ab) 
0.204 0.020 0.136 
0.480 0.048 0.320 
0.855 0.086 0.570 
1.460 0.146 0.975 
2.340 0.234 1.560 
2.295 0.295 1.965 


A, +A, 





TABLE 12 — NUMBER OF MOLECULAR LAYERS OF SODIUM HEXAMETAPHOSPHATE ON THE SURFACES 
OF 6 PER CENT M FRACTION WITH THE AXIAL RATIOS OF 1:100:200 


Number of Layers on the Surfaces of 





100 
Amount Added Amount Adsorbed 


Faces 
(Ib/bb!) Ib/bbI % A, (be) 


0.10 075 75.0 0.885 
0.25 0.183 73.0 2.080 
0.50 0.325 65.0 3.710 
1.00 0.556 55.6 6.320 
2.00 0.890 44.5 10.150 
4.00 1.120 28.0 12.800 





DECEMBER, 1961 


Az (ac) 


010 001 
Faces Faces 
A, (ab) 


0.443 0.004 0.295 
1.040 0.010 0.695 
1.855 0.019 1.235 
3. 160 0.032 2.110 
5.075 0.051 3.380 
6.400 0.064 4.270 


A, +Az 








Reduction in viscosity and gel strength was 
found to be directly related to the adsorption of 
polyphosphates. Adsorption of each polyphosphate 
indicates the same general mechanism of defloccula- 
tion. Most of the adsorption, (55 to 80 per cent) is 
believed to be chemical-bonding in nature and 
appears to be confined principally to the edge 
surfaces. Although some adsorption undoubtedly 
occurs on the flat surfaces, the evidence shows 
that such adsorption probably plays a minor role 
in effecting changes in the gel strength. 

Axial ratios of 1:15:20 are indicated for sodium- 
montmorillonite particles ranging from 0.1 to 0.04- 
micron average size. A 3:4 ratio of lateral dimen- 
sions appears correct. 
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Velocity-Log Interpretation: The Effect of 


Rock Bulk Compressibility 


J. GEERTSMA * 
MEMBER AIME 


ABSTRACT 


The relationship between porosity and the speed 
of propagation of acoustic waves in fluid-saturated 
porous rocks as measured by the Sonic log and by 
ultrasonic techniques is analyzed. Biot’s continuum 
theory lis used to explain the difference in acoustic 
wave propagation between a dry and a liquid- 
saturated porous material. The porosity is a vari- 
able in this theory. However, the acoustic wave 
propagation in the dry rock depends too on porosity, 
and this dependence is not predicted by the theory. 
Frequently in dry sandstones, a nearly linear rela- 
tionship between reciprocal acoustic wave velocity 
and porosity is observed in the low-porosity range. 
The physics behind this behavior is outlined. An 
empirical relationship of the form, 1/V =~ A +B ¢, 
applies accordingly for many porous dry rocks, 
provided the porosity is the only variable. The 
presence of a liquid in the pores changes the value 
of B, and this change is found to be in agreement 
with the Biot theory. 

The time-average relation introduced some years 
ago? results in an equation of the same type — 
1/V = d/V; + (1 — $)/V, — but is not based on a 
sound physical picture. Still, this relation some- 
times predicts approximately correct A and B values. 

Carbonate rocks with their complicated pore 
structures sometimes show a different relationship 
between wave velocity and porosity, unfavorable 
for log interpretation. Examples are presented. The 
simultaneous presence of calcite, dolomite and 
anbydrite, with their different grain densities and 
matrix compressibilities, complicates acoustic-log 
interpretation in carbonate rocks still further. 

Other complicating effects of acoustic-log inter- 
pretation are discussed. They are related to the 
influence of shale streaks and natural fractures on 
the average wave velocity observed by the logging 
tool and to the effect of adsorption phenomena on 
wave propagation in unstressed rocks particularly 
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in sandstones. 


INTRODUCTION 


The velocity of propagation of sound waves in 
porous sediments as a function of depth is a 
quantity frequently measured. Velocity loggers of 
various design may be used, for which velocity is 
defined as the distance between a wave generator 
and a wave detector (one-receiver system) or the 
distance between two wave detectors (two-receiver 
system), divided by the shortest time required for 
a vibrational pulse to cross this distance. Velocity 
loggers are used to assist in seismic prospecting, 
to differentiate between the various types of sedi- 
mentary rock layers, and also to determine rock 
porosity and pore fluid content. This paper is 
especially concerned with the relation between 
velocity and porosity and the effect of the pore 
fluid. 

Most velocity-log interpretation in terms of 
porosity is based on a formula advocated by Wyllie, 
et al,? suggested earlier by Hughes and Jones,* 
and known as the ‘‘time-average relation’’. This 
formula applies for a model consisting of a layered 
system of parallel alternating slices of a solid and 
a liquid, crossed by the wave path perpendicular 
to the interfaces. This must be an unsuitable model 
for the derivation of wave propagation properties 
of liquid-saturated porous media. It suggests that 
only rock matrix and fluid properties influence 
the wave velocity. The surprising fact, nevertheless, 
is that applications of this formula to clean sand- 
stones under specific conditions (sufficient effective 
stress) results in porosity predictions which are 
close to reality. 

Other authors, such as Gassmann,* White and 
Sengbush,5 Brandt,© and Hicks and Berry,” make 
use of more realistic models to arrive at a wave 
velocity formula for porous sediments. Their models 
consist of various liquid-saturated packings of 
frictionless spherical grains and are, therefore, 
also highly idealized. A number of experimental 
results obtained in the laboratory can be better 
explained with these model theories than by the 
time-average relation. A study of a general character, 
not depending a priori on a special model of the 
porous structure, was still lacking. 
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THE BIOT THEORY 


For the description of the speed of propagation 
of elastic waves in a two-component system, one 
has to consider two individual deformations. In a 
fluid-saturated porous system, these are the de- 
formation of the solid frame and that of the fluid. 
Following a two-dimensional study by Zwikker 
and Kosten8 on sound propagation in elastic porous 
media containing air, a three-dimensional general 
theory developed by Biot! showed that the de- 
termination of the relationship between the de- 
formations of solid and fluid is a problem which 
can be solved by the classical theory of me- 
chanics. 

If a porous solid is deformed, its pore volume in 
general also changes. If this pore volume contains 
fluid, this fluid can change in volume owing to its 
own compressibility. As the fluid compressibility 
in most cases will be different from that of the pore 
volume, the fluid will try to flow through the porous 
medium under the influence of an induced pressure 
gradient. A continuity equation applied to this 
case describes the balance between fluid displace- 
ment and change in fluid volume. The equations of 
motion of solid and fluid and the equations of state 
(the stress-strain relation of both components), 
together with the equation of continuity, are suffi- 
cient for the determination of the dynamic behavior 
of the porous system, provided the effect of temper- 
ature changes due to dissipation of energy as heat 
can be neglected. 

Applying this procedure to assumed isotropic and 
elastic porous media and further assuming that 
relative motion between rock grains and fluid takes 
place according to Darcy’s law, one arrives at a 
set of equations that describe the behavior of dila- 
tational waves (providing these waves are of a 
length greater than that of the largest grain dimen- 
sion). 


e2 


V 2 (He ~— KZ) = 
e ¢) 52 


(pe — p;¢), 


2 
V? (Ke -L¢) = SF (me ~— p.0)- a (1) 


In this set of simultaneous differential equations, 
e represents the absolute dilatation of the porous 
solid, and € the ‘“‘relative dilatation’ between 
solid and fluid. The proper definition of ¢ is 


C= ¢(e- 9, 


where ¢ represents the porosity and ¢ is the absolute 
dilatation of the fluid. Three densities appear in 
these equations: pis the total or bulk density of 
the system; p; is that of the pore liquid; and p, 
can be written in the form 
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where xk represents a mass-coupling factor which 
must have a numerical value > 1; x = 1 means that 
there is no mass-coupling. Zwikker and Kosten8 
obtained a theoretical value of x = 3 for a system 
of randomly oriented pores. We will consider x as 
an unknown material property. A method for the 
experimental evaluation of the mass-coupling 
factor will be published soon. y/k is the mobility 
of the pore fluid, H, Kand L are elastic deformation 
properties representing combined functions of the 
individual deformation properties of the components 
of the system and the porosity. 

The nature of the constants H, K and L can be 
analyzed with the use of a more detailed theory 
relating stress and strain in a fluid-saturated porous 
medium. In 1951, Gassmann” outlined the funda- 
mentals of this theory. Further studies on this 
subject have been published by Geertsma!° and by 
Biot and Willis11 Applying these published results 
to define H, K and L, we find (for derivation see 
Geertsma and Smit 12) 
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compressibility of rock matrix, 
cp, = compressibility of the rock bulk material, 
c,; = compressibility of the pore fluid, and 
G, = the shear modulus of the rock bulk material. 


The differential Eqs. 1 predict the existence of 
two plane dilatational waves with very different 
properties. Only the wave ‘‘of the first kind’’ is a 
true wave. The wave ‘‘of the second kind’’ obeys 
an equation which has more the character of a 
diffusion equation and is, thus, highly attenuated 
and only observable in the immediate neighborhood 
of the wave source.* For the interpretation of the 
velocity log, only the wave of the first kind has to 
be considered. The exact solution for the wave 
velocity of the wave of the first kind is easily 
arrived at but is impractical for engineering pur- 
poses. Geertsma and Smitl!2 derived a suitable 
approximate solution which reads 





*Patterson 13 reported rates of propagation which are very 
likely those of the wave of the second kind. Tchekaliuk !4 
reported some preliminary work on the theory of the wave 
of the second kind. 
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Vo represents the theoretical wave velocity at 
zero frequency, V,, is the wave velocity for infinitely 
large frequency and w is the frequency. Substitution 
of Eqs. 3, 4 and 5 into Eqs. 7 and 8 leads to 
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which, by taking an average value of Poisson’s 
ratio vz, = 1/5, reduces to 


8 (1 — B)2 ] 
y2 am B + —— a — . .(7b) 


And in the same way, 


B 4 
2 © al iby. 
wa-|(Es fo, 





Bd/ Kp, +(1-B)(1-B-2 d/«) l 
(1-¢-B)c, + dc, p} 
J pi-—— 
p «K 
tS Ee TE SS oe Seoe cae ee ee 


THE INFLUENCE OF FREQUENCY 


It is thus predicted that the wave velocity depends 
on frequency, and we have to know the extent of 
this dependence in practice. Since the differential 
Eqs. 1 are valid only so long as Darcy’s flow law 
is valid, the applicability of solution Eq. 6 is 
restricted to small values of w/w,. The high- 
frequency range has also been investigated by 
Biot, but it appears that in this range neither the 
value of V.. nor the general picture of the frequency 
dependence of wave velocity is altered. 





*The ratio @/®, , the frequency number, may be considered 
as the Reynolds number of the vibrating motion, It represents 
the ratio of inertia and viscous forces. 
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Eq. 6 predicts that for w/w, < 0.1 and for @/w, 
> 5 the wave velocity is nearly independent of 
frequency. In the range between these limits there 
is a notable dispersion. At first sight, therefore, 
one should conclude that, for an easy interpretation 
of wave velocity, the frequency of the signal should 
not be chosen in the ‘‘transition zone’’. In practice, 
however, this depends on the numerical differences 
which can be expected between Vp and V,,. 

The location of the transition zone is indicated 
in Fig. 1. Here, for the case of water-saturated 
sandstone, @ has been plotted against porosity 
with permeability as a parameter. Water viscosity 
was taken to be p = 1 cp, while for the densities of 
solid and liquid the values 2.5 and 1.0 gm/cc were 
taken, respectively. The formulas for the constant- 
permeability curves follow from Eq. 9. Assuming 
k = 3, one finds for the lower limit of the transition 
zone 


1-06 d\¢ 

ad po = 

low = 10 G << $) a 
while 


“high = 50 Mlow 


represents the upper limit of the transition zone. 
Only the curves for the lower limits are given in 
the figure for four values of the permeability. 
Though there is no theoretical relationship between 
porosity and permeability, we shall use a rough 
empirical correlation between the two quantities 
obtained from data collected by Muskat!> If this 
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FIG. 1— GRAPH INDICATING REGION 

WHERE WAVE VELOCITY IS FREQUENCY- 

DEPENDENT FOR WATER-SATURATED 

SANDSTONES (MASS-COUPLING FACTOR 
Kk = 3). 














correlation is applied, the frequency-dependent 
region can be roughly evaluated. This is represented 
by the L-shaped section in the figure. The frequency 
region used in practice in acoustic-logging tech- 
niques lies (see Striplingl!6) between 10 and 20 
kc/sec. This region is also indicated in Fig. 1. 
It is seen that the frequency-dependent region is 
not completely avoided, but for porosities below 
20 per cent the velocity approaches Vo, even for 
low liquid viscosities such as that of water. 

On the other hand, laboratory experiments with 
ultrasonic devices will measure the approximate 
value of V,,. A theoretical evaluation of the ratio 
V,. /Vo reveals that, when the mass-coupling factor 
kK approaches o, the frequency dependence dis- 
appears and V,, = Vp. In practice « is not infinite, 
yet a numerical study of the ratio V,,/Vo for sedi- 
mentary rocks reveals that the difference between 
Voo and Vo using reasonable values for « will in 
general be smaller than 10 per cent. A more precise 
determination of the differences between V,, and 
Vo can be given only after a correlation between 
B and ¢ has been discussed. In most cases the 
actual differences are not large, although they are 
not entirely negligible if one compares sonic and 
ultrasonic data. 

Fig. 1 also indicates the limit of validity of the 
continuum theory. This limit is given by the con- 
dition that the wave length should be much larger 
than the grain diameter. It is assumed here that 
for sandstones the limiting wave length is 0.5 cm. 
This condition leads (using data presented later 
on) to the limiting curves in the upper part of the 
figure. The conclusion can be drawn from these 
curves that the Biot theory is applicable up to 
frequencies far in the ultrasonic region. 


INFLUENCE OF THE BULK DEFORMATION 
PROPERTIES 


Eq. 7 will be used from now on for the interpre- 
tation of acoustic wave velocities in saturated 
sediments. It gives the effect of the properties of 
rock matrix, rock bulk material and pore fluid on 
the wave velocity. Solid and liquid densities, 
porosity and compressibilities of solid matrix and 
liquid occurring inthis equation can be considered 
as true and mutually independent constants. These 
same variables (though indirectly) appear in the 
time-average relation. In the present equation, 
however, the elastic properties of the rock bulk 
material (expressed by f) play a dominant role. 

The first term in the velocity equation in general 
provides by far the largest contribution to the 
velocity value and is a direct function of the rock 
bulk elastic-deformation properties c, and v,. An 
exception is found for suspensions, where one has 
to insert B = 0; accordingly, for this case the wave 
velocity depends only on porosity (or, rather, con- 
centration) and the densities andcompressibilities 
of the individual phases. 

There are sufficient indications that the rock 
bulk compressibility is a function of porosity, 
pore-size distribution, composition of the cementing 
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material, type of sedimentary rock and effective 
stress. 


DEFORMATION THEORIES BASED ON MODELS 
OF THE BULK MATERIAL 


The theory of Biot and its further development 
thus presents only a partial solution for the problem 
of acoustic-wave-velocity interpretation in terms 
of porosity. It predicts only the difference in 
acoustic wave velocity in dry and in fluid-saturated 
porous media. For further progress, basic knowledge 
on rock bulk compressibility is required. Many 
attempts to evaluate this property of sedimentary 
rocks on a theoretical basis can be found in the 
literature. Various models of porous structures 
have been analyzed for this purpose. Adams and 
Williamson!” as early as 1923 suggested two 
simple theoretical models for an investigation of 
the influence of pores on rock bulk compressibility 
in the low-porosity range. 

1. The unit element of the structure is a solid 
sphere with a central spherical cavity. 

2. The unit element is a solid sphere with a 
narrow spherical shell-shaped cavity, separating 
a solid central sphere from an outer spherical shell 
comprising the remaining solid material. 

We carried out the necessary calculations for 
these models using the classical theory of elasticity 
and found that, in the case of the spherical pores, 


M+@ 


“b = M( — ¢) °*? 


or 


G-)-FAA5- Lat gear aie ae 


M is only a function of Poisson’s ratio of the matrix 
material, 





2(1 — 2v) 
aie ae i 


This result means that (1/8 — 1) is, roughly, directly 
proportional to ¢ in the low-porosity range. It can 
be shown that, for any model of a solid with suffi- 
ciently isolated holes distributed at random, a 
relationship of the form 


ee. (10a) 


holds. Using van der Knaap’s!® definition of 
pore compressibility 


Cp = (1/B = 1) g/d, 


this implies that 


es lin ele re (11) 
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or a nearly constant pore compressibility in the 
low-porosity range. 

Conversely, for the model with the shell-shaped 
pores (shell crack at one-half the outer radius of 
the unit sphere), one obtains 


8M + 1 
Ch=—>TH Cr? 


‘G- ; (M+), 


So here, (5- ) is independent of porosity. This 





may possibly be a suitable model for a rock con- 
taining microcracks. A model whose deformation 
characteristics agree in some important aspects 
with those of granular sediments is a packing of 
spheres. Several types of packing are possible, 
and different quantitative results are consequently 
obtained by various authors. Apart from this de- 
pendence of deformation characteristics on packing, 
there is the problem that an exact calculation of 
the bulk deformation properties of packings of 
spheres is complicated. Most investigators avoid 
some of the complications by assuming that the 
spheres of the model deform without friction at 
the contact points. This means that at the grain 
contacts only normal forces are taken into account. 
Only Mindlin19 has considered the additional 
influence of friction. 


Two characteristic deformation properties of 
packings of spheres are the following. 

1. The elastic ‘‘constants’’ of an array of spheres 
are pressure-dependent. Calculations show that the 
bulk compressibility of such models is invariably 
proportional to (@ — p)71/3, 

2. The porosity of packings of uniform spheres 
varies roughly between 25 and 45 per cent, but 
from deformation studies of different uniform 
packings one cannot obtain a clear picture of how 
the compressibility of a uniform sphere packing 
depends on porosity. Brandt® did some interesting 
work on packings built up from spheres of various 
diameters. In this way the porosity could be varied 
more-or-less gradually. He found that the pore 
compressibility of nonuniform sphere packings 
should be nearly independent of porosity, as in 
the case of a solid with isolated holes. 


We conclude from this analysis that for many 
porous sediments, other variables kept constant, 
a relationship of the form of Eq. 11 possibly exists 
between pore compressibility and porosity. 


Introducing a constant C in Eq. 11, we can write 


$2 
a oi 2 « wae 


Ch = Cp + PCy = Cy + 


The dilatational wave velocity in the dry rock 
depends on this rock bulk compressibility. 


DECEMBER, 1961 








2 {(%) {(%) 


— €BPr (1-6) og -#)|e, + 13) 


Assuming as a first approximation that Poisson’s 
ratio is independent of porosity, we realize that 
K(vb)/prce = Ves the matrix wave velocity. So 


1/v = 1/V, J (ES -1)¢ 


which for small porosity values can be further 
approximated as 


1/V = 1/V, frt(- ) ‘| ; 


This equation is of the form 








ee ee 


and the additional information is that A = 1/V, in 
this approximation. The coefficient B depends 
strongly on the value of the pore compressibility, 
but it has nothing to do with the presence of a 
fluid in the pores. Substitution of Eq. 12 in Eq. 7a 
for the case of a liquid in the pores does not 
change the approximate result Eq. 13 very much. 
To a first approximation, only A and B obtain a 
different value. 


EXPERIMENTAL DATA 
CLEAN SANDSTONES 


In Fig. 2, rock compressibility data calculated 
from velocity measurements for a number of dry 
sandstones at atmospheric pressure have been 


collected. They are derived from measurements 
by Wyllie, et al,2 data collected by Wuerker2% 
and a number of experiments carried out in our 
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FIG. 2 —c,/c, AS A FUNCTION OF POROSITY FOR 
DRY SANDSTONES ACCORDING TO DYNAMIC MEAS- 
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laboratory. In this figure, 1/B = (c,/c,) is plotted 
against the porosity. Results roughly suggest a 
nearly linear relationship 


(7/8 ~ 1) @ SOO. ks ck ses . .(10b) 


A relationship of the form of Eq. 10a would be 
more satisfactory from a theoretical point of view, 
but would be less convenient to handle while it 
hardly gives a better result for low porosities. 
Thus, for these sandstones, the pore compressibility 
is indeed nearly independent of porosity, in agree- 
ment with the foregoing theoretical considerations. 

For the dependence of f on effective stress, we 
can use van der Knaap’s18 empirical relationship 
for the bulk compressibility of sandstones, based 
upon an observed pore compressibility 


Cp = te) 3.5x10-3 (6 p)-2/3sq m/MN.* 


If it is assumed that c, = 0.25 x 1074 sqm/MN, 
this leads to 


i ) = 140¢ (79—p)-2/3. «2... .(14) 


Eq. 14 cannot be the general relationship, since 
it predicts an infinitely large bulk compressibility 
at (@ — p) = 0. This is only true for loose sand. A 
closer look at van der Knaap’s results, which are 
reproduced here in Fig. 3, shows that the straight 
line between log cy and log (@ — p) breaks down 


“4 sq m/MN (square meters per meganewton) = 6.89 x 10°3 
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FIG. 3 — PORE COMPRESSIBILITY DATA OF SAND- 
STONES (ACCORDING TO VAN DER KNAAP IN REF. 
18); 1 MEGANEWTON/SQ M = 145.04 PSI. 








somewhere below (@ — p) = 2 or 3MN/sq m. At low 
effective stress values, the curve bends down to 
a horizontal asymptote; the pore compressibility 
for (© — p) = 0. Additional experimental work has 
confirmed this, and in Fig. 4 the final result is 
indicated. This result suggests that for consolidated 
clean sandstones an empirical relationship of the 
form 

a 

n 


= ee 2 : ~ 
G i * fo (ap) + (ag) ] *~ -as) 


will be the best average expression for the bulk 
compressibility. If we take a = 50, n = 3/2 and 
b = 5,000 (MN/sq m)”%, one obtains Eq. 10b for 
(G — p) = 0, while Eq. 14 is approached for large 
values of (@ — p), since 


5,000-2/3 = 3.5 x 10-3, 


The coefficients a, 6 and n will depend on grain 
cementing material, pore configuration, etc., but 
at any given stress level there is a linear relation 
between the porosity and (1/8 — 1). 

Combination of Eq. 15 with the general theory 
just outlined permits a more precise determination 
of the ratio V,,/Vp in sandstones. This ratio is 
plotted in Fig. 5 for an average value of the mass- 
coupling factor x = 3 and for water-saturated sand- 
stones at different values of (@ — p) and @. It is 
seen that the ratio V,, /Vp is in general lower than 
1.0. 

It is now possible to calculate the velocity of 
plane dilatational waves as a function of porosity 
at zero effective stress in water-saturated sand- 
stones according to Eq. 8a. The results are given 
in Fig. 6 where V,, has been plotted for x = 3. 
These theoretical results are compared with the 
ultrasonic measurements reported by Wyllie, et al, 
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and it is seen that the agreement is good over the 
whole porosity range. 

Next, V,, and Vp in water-saturated sandstones 
can be calculated as a function of both porosity 
and effective stress. Results are plotted in Fig. 7. 
In this figure the time-average relation is also 
indicated. A single curve is obtained for the time 
average relation since it is independent of effective 
stress. It is seen that, although the difference 
between the results of the two theories is con- 
siderable for moderate porosity values at very low 
and high effective stresses, for the stresses pre- 
vailing at average reservoir depths the time- 
average relation may predict more-or-less accurate 
porosity values. 

Eq. 15 also explains the entirely different velocity 
data obtained by Hicks and Berry7 and by Wyllie, 
et al.2 The measurements by Hicks and Berry were 
cafried out in water-saturated sandstones at an 
effective stress of about 42MN/sq m. In Wyllie’s 
measurements, the effective stress was zero. Fig.8 
demonstrates that velocity predictions with Eq. 15 
for (9 — p) = 50 MN/sq m agree very well with 
Hicks and Berry’s results, while those predicted 
for (G — p) = O fit in quite well with Wyllie’s data 
(as was already shown in Fig. 6). Fig. 6 shows 
that the presence of water in the pores straightens 
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FIG. 5 — THE RATIO Y,/V) FOR WATER-SAT URATED 
SANDSTONES USING THE EMPIRICAL RELATION, 
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FIG. 6 — COMPARISON BETWEEN MEASURED AND 
CALCULATED ULTRASONIC PLANE DILATATIONAL 
WAVE VELOCITIES IN AIR-DRY AND WATER-SATU- 
RATED SANDSTONES (MEASURED DATA ACCORDING 
TO WYLLIE, ET AL). 


the curves of 1/V against ¢, and an approximate 
relationship of the form 


1/V=A+B¢ 


applies over a larger porosity range in the water- 
saturated situation. More detailed information on 
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TABLE 1 — SOME AVERAGE VALUES FOR MATRIX 
DENSITIES, MATRIX COMPRESSIBILITIES AND WAVE 
VELOCITIES FOR CALCITE, DOLOMITE AND ANHYDRITE 


Plane 
Dilatational 








G:ain Matrix Acoustic Wave 
De.isity Compressibility Velocity 
Pr c, x 1074 m?/MN (m/sec) 
Calcite 2.70 —2.76 0.135 6400 
Dolomite 2.81—2.86 0.100 7200 
Anhydrite 2.96—3.00 0.135 6000 





the difference between the acoustic wave velocity 
in dry and in water-saturated sandstones is given 
in Fig. 9. Here, the ratio y—" ig plotted 
against the porosity, in which it is assumed that 
(1/B — 1) = C¢ and C is taken as a variable. Low 
C-values apply to tight rocks or rocks under high 
confining pressure. For Vg the case of water in the 
pores is considered (p; = 1.00 and c; = 4.00 x 
10-4 sq m/MN). It is seen that for low C-values the 
effect of water in the pores becomes small and 
that the difference Vo — Vgyy can even be negative, 
that is, the water in the pores lowers the wave 
velocity. It is concluded that acoustic velocity 
data are able to differentiate between gas or liquid 
in the pores in soft formations only. 


CARBONATE ROCKS 


The different origin of carbonate porosity, to- 
gether with the mixed presence of calcite and 
dolomite sometimes contaminated with anhydrite 
inclusions, suggests a priori that acoustic-log 
interpretation in carbonate rocks in terms of porosity 
will be more complicated than in sandstones. 

Calcite, dolomite and anhydrite have different 
matrix densities and matrix compressibilities. Some 
average values are collected_in Table 1. 

Therefore, from the beginning we will make a sub- 
division in calcite and dolomite and exclude the 
(in general, nonporous) anhydrites. We determine 
ultrasonic wave velocities in a series of carbonate 
rock samples of which van der Knaap!8 has reported 
the pore compressibilities. We will refer to this 
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FIG. 9 — EFFECT OF THE PRESENCE OF WATER ON 
ACOUSTIC WAVE VELOCITY ACCORDING TO BIOT’S 
THEORY. 
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series as Series A. In Ref. 18 van der Knaap showed 
that for these carbonates the pore compressibility 
was strongly influenced by the porosity. At a given 
stress level, cy increases with decreasing porosity, 
particularly in the low-porosity range. If our theo- 
retical conditions are correct, this predicts a quite 
different relationship between acoustic wave veloc- 
ity and porosity than observed for sandstones, at 
least in this low-porosity range. Because van der 
Knaap’s core material was partly calcite and partly 
dolomite, the velocity data had to be subdivided 
necessarily according to the type of carbonate. 
Unfortunately, this also resulted in a subdivision 
of the porosity range — the calcites were predom- 
inantly low porous, the dolomites showed higher 
porosity figures (9 to 17 per cent). In agreement 
with the pore compressibility data, the wave veloc- 
ities varied considerably with effective stress. 
In Fig. 10, therefore, reciprocal wave velocities in 
the calcites are plotted as a function of porosity 
for the arbitrary stress level of 50 MN/sq m (about 
500 atm). In Fig. 11 the corresponding data for 
the dolomites are shown. In neither plot can a 
linear relationship between 1/V and ¢ be detected; 
the reciprocal velocity varies rapidly in the low- 
porosity range and slowly at higher porosity. This 
is exactly what should be expected on the basis of 
the pore compressibility behavior. 

Not all carbonates are alike, however. Also 
plotted in Figs. 10 and 11 are data obtained from 
dry cores of another carbonate reservoir, determined 
at the same effective stress of SOMN/sq m. Again 
calcites are separated from the dolomites. This 
series, indicated by B, does show a nearly linear 
relationship and thus behaves similarly to many 
sandstones. A comparison between Curves A and 
B in the same figures further reveals that wave 
velocities in Series A for a given porosity are 
invariably lower. Accordingly, Formation A has 
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to show the largest pore compressibilities. Addi- 
tional pore compressibilities measured by van dex 
Knaap for Series B are reproduced in Fig. 12. Com- 
parison between these data and the data pertaining 
to Series A reproduced in van der Knaap’s publica- 
tion!8 reveals that, indeed, Series A has by far 
the larger pore compressibilities, particularly in 
the low-porosity range. It is also apparent from 
Fig. 12 that in Series B the pore compressibility 
is independent of porosity, in agreement with the 
linear relationship between 1/V and o. 

In Fig. 13, finally, acoustic velocity data of 
various calcites are compared. Sonic-log data 
obtained in a typical calcite section in an Austrian 
well at an average depth of 1,400 m are compared 
with ultrasonic wave velocities in the calcite Series 
B just mentioned, at a corresponding effective 
stress of 150 atm. The agreement is satisfactory. 
In addition, some ultrasonic velocity measurements 
in West European calcite-outcrops are included. 
These data also seem to agree with the other data. 
For the sake of completeness, the time-average 
relation assuming the presence of water in the 
pores is also given. This line fits the data as well. 
It must be noted, however, that the laboratory 
measurements were taken on dry cores. Again we 
are forced to conclude that the time-average relation, 
although untenable from a physical point of view, 
can by a lucky coincidence give correct results. 


INFLUENCE OF SHALE LAYERS AND 
HORIZONTAL FRACTURES ON EFFECTIVE 
VERTICAL STRESS 


We have seen that a decrease in bulk compress- 
ibility with increasing effective stress is a property 
of most sedimentary rocks. It causes an over-all 
gradual increase in elastic wave velocity with 
depth. At the same time, any stress discontinuity 
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FIG, 11 — ULTRASONIC WAVE VELOCITIES IN DOL- 
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in the neighborhood of the borehole must affect 
wave velocity. Such discontinuous wave velocities 
might easily be interpreted as porosity variations 
if one were not aware of the effective-stress effect. 
Three examples are given as follows. 

1. So long as the fluid contact in the pores is con- 
tinuous in a downward direction, the gradient of 
the effective vertical stress is fairly low. Imperme- 
able shale layers, however, may sustain a consider- 
able fluid pressure difference at both sides of the 
permeability barrier. Even if this impermeable 
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barrier separates two parts of the same porous 
sediment, the wave velocity in these parts may 
be different owing to different effective-stress 
conditions. 

2. Shale layers may disturb the vertical stress 
distribution in yet another way. Zheltov and Khris- 
tianovich,2! in a detailed study, showed that a 
considerable re-distribution of vertical effective 
stress may be caused by plastic clay formations 
of little strength. A plastic stress distribution 
around the well will result in a stress relief of 
the sediments immediately below the clay layer, 
while the rocks lying above the clay form a dome 
and experience an increased effective stress. If 
the layer below the clay is a producing formation, 
the low wave velocity found at the top of this 
formation might easily be interpreted as a high 
porosity. That such plastic stress conditions in 
clays around the well are real has been shown by 
Hicks.22 Using a multi-receiver velocity-logging 
system, he found that the wave velocity in a 
clay layer increased with increasing receiver- 
transmitter spacing. The larger this distance, the 
deeper the penetration of the sound waves. This 
showed that there was a considerable pressure 
relief in the clay near the borehole. A further con- 
sequence, not mentioned by Hicks, is that this 
stress relief makes itself felt immediately below 
the shale layer. The distance over which it is 
felt outside the clay depends on the extent of the 
plastic zone in the clay, and this in turn depends 
on the deformation properties and strength charac- 
teristics of the clay and on the depth of the layer. 

3. The same effect, but still more pronounced, 
may occur if a horizontal or nearly horizontal 
natural fracture intersects the wellbore. In the 
fracture, only the fluid pressure is available to 
sustain the vertical load of the overburden. This 
causes stress concentrations around the fracture, 
the formation directly above the fracture experiences 
an increased effective stress and the formation 
immediately below the fracture is relieved. 

The resulting stress pattern leads to a reduction 
in velocity in the formation below the fracture 
which is greater than the increase in velocity in 
the formation above the fracture. The average 
result, ‘‘seen’’ by the log is a velocity reduction. 





PHYSICAL INTERACTION OF PORE FLUID 
AND POROUS SKELETON 


van der Knaap!8 and others have noted that the 
nature of the fluid, and also the amount of fluid 
present in the pores, may have a remarkable effect 
on the elastic properties of sedimentary rocks, 
This behavior results in experimental data on 
samples, particularly at atmospheric pressure, 
that are completely different from the theoretically 
predicted results, and one can easily misinterpret 
the observed wave velocities. 

This will be illustrated with some examples. 
The resonance frequency (being directly proportional 
to wave velocity) of porous rods was measured for 
both torsional and longitudinal vibrations. The 
cylindrical samples were examined both air-dry and 
saturated with different liquids at atmospheric 
pressure and under frame pressure. In the latter 
case, the sample was covered by a very thin rubber 
jacket, placed in a pressure vessel and pressurized 
on the outside by means of hydrogen. 

At atmospheric pressure the resonance frequen- 
cies of sintered porous material show a small 
decrease upon saturation with fluid. This is in 
agreement with the theory previously outlined. 
By contrast, corresponding measurements on sand- 
stones show quite different results in this respect. 
This is illustrated by the data collected in Table 
2. The figures in this table give reasonance fre- 
quencies for a number of sandstones and sintered 
porous structures saturated with various liquids. 


In the case of sandstones, only the addition of 
heptane and sym-tetrachloroethane results in a 


small decrease in resonance frequency. For sintered 
materials, all liquids show this effect. However, 
saturating the sandstone samples with water clearly 
results in a much greater drop in resonance fre- 
quency. Saturating the samples with oils gave rise 
to a considerable increase in resonance frequency, 
compared with that of the dry samples. One easily 
arrives at the conclusion that the wave velocity 
should make a clear distinction between oil-bear- 
ing and water-bearing formations. These facts 
cannot be explained by the Biot theory. The exper- 
iments also showed that the attenuation of elastic 
waves in water- and oil-saturated sandstones, as 
derived from the resonance tests, is considerably 





TABLE 2 — RESONANCE FREQUENCIES OF WET AND DRY POROUS SYSTEMS AT ATMOSPHERIC PRESSURE 


















esas Resonance Frequency, Wet System 
Dry Sym-tetra- Transformer Middle East 
k System * Heptane chloroethane Oil Crude Water 
Porous System (per cent) (darcies) 1 t 1 t 1 t 1 t 1 t 1 t 
Rumaila Ss 20 ~ 12900 8600 12550 8450 11250 9000 13350 8600 14450 8950 9180 7000 
Zubair Ss 17 - 11600 8750 11200 8550 11550 8700 13550 8950 14100 9550 7600 5600 
Berea Ss 24 0.24 12150 8750 11450: 8300 11200 8200 14000 9400 14300 9075 7400 5800 


Sintered Glass 36 39.5 10300 7150 9900 7200 


Si Carbide 26 1.9 22400 14600 21700 14100 
Carborundum 31 14.3 31700 20500 31000 20100 
Heptane : p = 0.68, pp = 0.40 cp 


Sym-tetrachloroethane: p = 1.50, pp = 0.76 cp 


* 1= longitudinal waves, t = torsional waves. 
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- - 9600 6650 - - 9850 6750 
_ - 21600 14000 - - 21400 13950 
30400 19500 30800 19800 ~ ~ 30000 18500 
Transformer Oil pL = 25 ep 
Middle East Crude : pp = 71.5 cp 
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higher than predicted by the Biot theory. This 
derivation is omitted here but can be found in a 
paper by Geertsma and Smit. !2 

The observed anomalous behavior very probably 
is due to adsorption of the pore fluids and seems 
to be most important in an analysis of the physical 
behavior. of the so-called weathered layer. The 
reason why there is no need to explain these experi- 
mental results in the present study is that these 
adsorption effects do not show up in the practice 
of velocity logging. This was concluded from a 
series of experiments which showed that these 
anomalous effects become almost unnoticeable 
and that results are always in line with what should 
be expected from Biot’s theory provided a suffi- 
ciently large effective stress is applied to the 
skeleton.* 

Except in the shallow depth range, the effective 
stress around the borehole is sufficiently high to 
eliminate adsorption effects. We were careful not 
to include data influenced by such adsorption 
phenomena in the present study. The interesting 
effects of adsorption and its influence on the 
mechanical properties of porous rocks are still 
subjects of investigation. 


CONCLUSIONS 


This investigation of the relationship between 
the speed of propagation of elastic waves in liquid- 
saturated porous rocks and porosity leads to a 
number of conclusions. 

1. Unfortunately, the time-average relation now 
in use for deriving porosity values from the acoustic 
log is based on a wrong physical picture. However, 
in many sandstones and carbonates it predicts, 
independently of the fluid content of the pore space, 
about the correct trend of the relationship between 
acoustic travel time and porosity at depths between 
1,000 and 3,000 m, if for the wave velocity of the 
liquid in that equation the velocity in water is 
taken. 

2. The Biot theory explains the difference in 
acoustic wave velocity in dry and liquid-saturated 
porous mediaunder reservoir conditions of effective 
stress. This difference is usually rather small in 
deeply buried well-cemented reservoir rocks. Dif- 
ferentiation between liquids and gas with the Sonic 
log is possible only in near-surface rocks, preferably 
consisting of loose particles. 

3. In the dry rock frequently a nearly linear 
relationship between reciprocal acoustic wave 
velocity and porosity is observed, corresponding 
with a pore compressibility which is independent 
of porosity. The basic relationship between the 
deformation properties of the rock bulk material 
and the acoustic wave velocity is shown to be 
equal to that in nonporous solids. 

4. Carbonate rocks sometimes show a rather 
peculiar relationship between pore compressibility 
and porosity, as has been reported earlier. This 





*Adsorption effects also disappear if the frequency is chosen 
in the higher ultrasonic ‘region above about 1 megac ycle/ 
second, 
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corresponds with an equally peculiar relationship 
between acoustic wave velocity and porosity. 

5. The interpretation of the acoustic log in 
carbonate rocks is furthermore complicated by the 
presence of both calcite and dolomite, which have 
different matrix velocities. 

6. Both compressibilities and acoustic wave 
velocities are nearly always strongly pressure- 
dependent, Therefore, porosity-velocity correlations 
must always be carried out at corresponding stress 
conditions. 

7. The acoustic wave velocity in liquid-saturated 
porous sediments is somewhat frequency-dependent. 
The effect of frequency results in minor numerical 
differences (less than 5 per cent) between sonic- 
and ultrasonic-wave velocity data. 

8. The fact that acoustic wave velocity is a 
function of effective stress renders the Sonic log 
sensitive to discontinuities in the vertical effective 
stress distribution around the borehole. Such stress 
discontinuities may result from impermeable barriers 
between permeable layers and plastic stress con- 
ditions in clay layers around the borehole; they 
also may be caused by natural fractures. Velocity 
variations resulting from these phenomena can 
easily be misinterpreted as being due to porosity 
variations. 

9. The results of wave propagation experiments 
in sedimentary rocks under low or zero effective 
stress often are influenced by adsorption of the 
pore fluid on the grain surfaces. In that case, they 
are of no value for the interpretation of what the 
Sonic log ‘‘sees’’ at depth. 


NOMENCLATURE 
c, = rock bulk compressibility, m=! L ¢2 


c, = liquid (fluid) compressibility, m7! L t? 
c, = pore compressibility, m= ?L t2 
c, = rock matrix compressibility, m=! L t2 
e = absolute dilatation of rock bulk material 
G, = shear modulus of rock bulk material = 
3(1_ — 2v) of, s 
2c,(1 + v) ’ mL e~ 
H,K,L = deformation constants occurring in Biot’s 
theory, m L-1 ¢-2 
k = permeability, L? 
p = hydrostatic pore tension, m L-! t~2 
t = time, t 
V = acoustic wave velocity (real part), L t~! 
Vo = acoustic wave velocity for zero frequency, 


Le 
V., = acoustic wave velocity for infinite high 
frequency, L t7} 
B=G/cp 


~ 
i 


absolute dilatation of pore fluid 


ww 
" 


relative dilatation between solid and 
fluid = d(e-«) 
K = mass-coupling factor > 1 


p = fluid viscosity, m L-! ¢-} 












v, = Poisson’s ratio of rock bulk material 
p = total or bulk density, m L-3 
Pp, = fluid density, m L-3 
Pp . 
Peo = ee m L 3 
@ = external hydrostatic pressure, m L-1 ¢ -2 
= porosity 
@ = angular frequency, t~! 
@®, = characteristic angular frequency = 
: -— 
md 5 tel 
RP PPc — Pi 
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DISCUSSION 


Vv. S. TUMAN 
MEMBER AIME 


Geertsma’s analytical work is very interesting, 
but we wonder if the idealized average stress 0 
that he has used in his study can be justified for 
the wellbore conditions. 

As he has mentioned, there will be stress relief 
at the wellbore, but these will appear in the shales 
as well as in the porous formations. Abnormal or 
residual stresses, not related to the overburden 
pressure, exists in the folded formations and in 
the vicinity of salt plugs. The temperature gradient 
between the wellbore and formation, due to mud 
circulation, will impose another set of stresses 
not related to the overburden pressure. 
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In permeable beds, a hydrostatic pressure gradient 
of several hundred pounds will also exist between 
the mud column and the formation fluid. Therefore, 
a ‘variable effective stress’? exists, which in- 
creases radially from the wellbore. 

In such circumstances, due to the stress envelope, 
sonic beams entering the formation with angles 
smaller than the critical angle are continuously 
refracted; following a curved path, they finally 
arrive back at the wellbore. 

Such beams with a curved path, with some specific 
assumptions, arrive at the wellbore before the 
refracted beams with critical angle. The beams 
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with the curved path are affected by the vertical 
and horizontal effective stresses. Since the 
horizontal stresses differ considerably from the 
vertical stresses, we wonder if an average stress 
ie GQ +p + © 
— sh 





can be applied under the wellbore 


conditions. 

According to Topping,* for a semi-infinite elastic 
solid where lateral deformations have been pre- 
vented, the radial stress relief is given as 


a2 a a2 
= — Br * Pe (75) (4-)) 


a = radius of the well, 
py = mud density, 
p = bulk density of formation, 
z = depth of the formation, 
v = Poisson’s ratio, and 
r = radial distance from the wellbore. 
Notice that at the equilibrium condition the hor- 


oe 2 
izontal stress is given by a, =3~ when v = 0.25, or 





* Topping, A. D. and Miles, A. J.: ‘Stresses Around a Deep 
Well’’, Trans., AIME (1949) Vol, 179, 186. 


AUTHOR’S REPLY 


Under field conditions in not too shallow inter- 
vals, laboratory tests indicate that the speed of 
propagation of plane dilatational elastic waves is 
indeed a function of the effective stress in the 
direction of propagation only. This is another con- 
sequence of the nonlinear elastic behavior of porous 
sedimentary rocks. These rocks, even when isotropic 
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the radial stress in one-third of the vertical stress. 

The velocity of sonic propagation is probably 
a function of the effective stress in the direction 
of propagation. This we believe is illustrated by 
two experiments performed by Wyllie, et al (Figs. 
2 and 3).** 

Wyllie and co-workers observed that, when a cirt- 
cumferential pressure was applied to the sample, 
the velocity of sonic propagation along the axis 
was less than when similar pressure was applied 
hydrostatically, or along the axis only. However, 
Wyllie offers a different explanation. 

Since the Topping relation is purely theoretical, 
some further experiments along the lines of those 
performed by Wyllie and his colleagues would be 
extremely useful. 

With this question in mind, we wonder if one 
should solve the differential equation of the wave 
propagation in terms of bulk compressibility, or 
whether it should be confined to the effective 
strains in the direction of propagation. We are 
inclined to think that the latter choice is more 
desirable. 





** Wyllie, M. R. J., et al: ‘An Experimental Investigation of 
Factors Affecting Elastic Wave Velocities in Porous Media”’, 
Geophysics (1958) Vol. 23, 459. 


TO V. 8S. TUMAN 


in the unstressed state, become anisotropic under 
nonuniform stress conditions. 


Some years ago, Wyllie, et al,** reported some 
comparative measurements of ultrasonic elastic- 
wave velocities under various stress conditions, 
i.e., uniform, axial and circumferential pressure 
exerted on a jacketed cylindrical sandstone sample. 
Such measurements are required for translating 
laboratory experiments under, say, uniform pressure 
into data pertaining to field conditions. Wyllie 
reported results for extreme stress conditions, the 





























Ultrasonic wove velocity, V 
$000 —— 
} 
Oox 
er 500 
4800 a ro ie 
Saxe Sra jn 300 
-- 
ee, 
ax. ——¢ 200 
a - 
i - 
a 
| »Z. 2inn 
4600} rai 1 ina 
| 4 
/ 
| / b id > w' —_ 
wot _| | 
| | 
| ——_ De constant 
——— Oa Srog 
| | 
42004 l i j 
C) 100 200 300 500 _ 
Radia! stress, 


FIG. D-2— VARIATIONS IN ULTRASONIC WAVE 
VELOCITY UNDER TRIAXIAL STRESS CONDITIONS 
— SAHARA SANDSTONE, POROSITY 8 PER CENT. 


247 





case of axial stress without circumferential stress 
and the effect of circumferential stress without any 
axial stress. 

To study the effect of nonuniform stress on elastic 
wave velocity in more detail, we performed a series 
of experiments with jacketed samples in a triaxial 
loading apparatus. Both axial and circumferential 
pressures could be varied systematically. The wave 
velocity was always measured in the axial direction. 
Results obtained on two sandstone samples — a 
Mainzer outcrop sample of 18 per cent porosity and 
a Sahara reservoir core sample of 8 per cent porosity 
— are shown in Figs. D-1 and D-2, respectively. 
Families of curves showing ultrasonic wave velocity 
vs radial stress are shown for various axial stresses. 
The curves for uniform pressure (0g, = Oqq), con- 
structed by connecting points of equal axial and cir- 
cumferential pressure are also shown. Raising the 
circumferential pressure above the axial pressure 
does not harm the samples. If the circumferential 


pressure is reduced considerably below the axial 
pressure, however, the sample may be crushed, and 
we stayed within safe limits in exerting this latter 
type of loading. Our results are more comprehensive 
than the curves reported by Wyllie, et al. Above 
100 atm, the radial stress has no influence at all] 
on the wave velocity. This follows from the fact 
that all experimental points for o,, 2 100 atm are 
of horizontal lines. 

For lower axial loadings, the circumferential 
stress influences the axial wave velocity to some 
extent. This effect is larger for the Mainzer sand- 
stone than for the Sahara core. In practice, these 
results show that, as far as the acoustic log is 
concerned, it is primarily the vertical stress that 
governs the velocity and not the mean hydrostatic 
stress. Therefore, we agree with Tuman that in a 
more refined theory we must take these anisotropic 
effects into account. However, the present paper 


was meant to outline only the basic principles. 
kkk 
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Additional Thermal Data for Porous Rocks--Thermal 


Expansion and Heat of Reaction 


WwW. H. SOMERTON 
MEMBER AIME 
M. A. SELIM 


ABSTRACT 


Thermal expansions and heats of reaction of 
three typical sandstones were measured in the 
temperature range of 25° to 1,000°C. The significance 
of these data in subsurface heat-transfer calculations 
is discussed. 

Linear thermal expansions were measured both 
parallel and perpendicular to the bedding planes. 
Volume expansions are reported as the perpendicular 
expansion plus two times the parallel expansion. 
Expansion behavior of the sandstones was found 
to be controlled by the expansion characteristics 
of the quartz content. Differential expansion of 
the quartz grains and other minerals included in 
the sandstones caused permanent deformation of 
the heated samples after they were cooled to room 
temperature. Structural damage resulting from heat- 
ing is probably an important cause of the reduction 
of thermal conductivity of heated samples. 

Measurements of heats of reaction were based 
on the known heat required for a-[3 quartz inversion. 
Thermal reactions, which probably include deby- 
droxylation of clay minerals and decomposition 
of carbonate minerals contained within the samples, 
were found to require more than one-fourth of the 
amount of heat necessary to raise the temperature 
of the rock alone. In shales and limestones, the 
reaction heat could be substantially greater than 
that required from specific heat considerations 
alone. 


INTRODUCTION 


In earlier work,1-3 methods of measuring thermal 
properties were developed and thermal data for 
several sedimentary-rock types were reported. Data 
included thermal conductivity, specific heat and 
thermal diffusivity. In addition, calculations of 
heat transfer in subsurface formations may require 
data on the thermal-expansion characteristics of 
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rock and the heats of reaction of mineral constit- 
uents. 

Thermal expansion of rocks is relatively small 
in magnitude and, from the standpoint of change 
in bulk density, has only minor effects on heat- 
transfer characteristics. Thermal-expansion be- 
havior, however, may have significant effects on 
the structure of rocks. Differential expansion of 
different minerals and along different crystal- 
lographic axes may result in structural damage which 
could effectively reduce the thermal conductivity 
of the rock. 

Many mineral constituents of sedimentary rocks 
undergo phase or related changes when heated to 
sufficiently high temperatures. Although some of 
these reactions (such as the a-f quartz inversion) 
are reversible and the heat absorbed is returned to 
the system upon cooling, many of the major reactions 
are irreversible. Through certain temperative ranges, 
the additional heat energy required to complete the 
reactions may be nearly as great as that necessary 
to raise the temperature of the rock in the absence 
of thermal reactions. If the heat requirements for 
these reactions are not considered, serious errors 
in heat-transfer calculations may result. 

Thermal-expansion and thermal-reaction data 
have been obtained for three typical sandstones. 
Methods of measurement and results obtained are 
reported and discussed in the following. 


THERMAL EXPANSION OF SANDSTONES 


The linear expansions of three outcrop sandstones 
(Bandera, Berea and Boise) have been measured 
in the temperature range of 25° to 1,000° C. Expan- 
sion measurements were made in the directions 
parallel and perpendicular to the bedding. 

The differential expansion apparatus used is 
that described by Mitoff and Pask.* The test sample 
is heated in an electric furnace with a temperature 
rise of 6°C/minute. The lengthening of the samples 
upon heating is compared with the small and known 
expansion of a fused silica rod. The change in length 
is transmitted to an X-Y recorder by means of a 
Statham transducer. A maximum error of + 1.5 
per cent has been obtained with this apparatus 
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for materials of known linear expansion. Upon 
reaching maximum temperature, the sample is 
cooled at approximately the same rate and the con- 
traction is recorded. The final length of the sample 
is measured to test the reliability of the final 
recorded value. 

The linear expansions of Berea sandstone both 
parallel and perpendicular to the bedding are 
compared in Fig. 1 with the expansion of quartz 
perpendicular and parallel to C-axis.> The most 
important features of Fig. 1 are: (1) the close 
agreement between the heating curves for Berea 
sandstone and the curve for quartz perpendicular 
to the C-axis; and (2) the lack of agreement of the 
Berea heating curves with the curve for quartz 
parallel to the C-axis. Discontinuities in the curves 
at 575°C are attributed to the a-f inversion of 
quartz. The cooling curves for Berea deviate con- 
siderably from the heating curves, the samples 
being permanently elongated at the conclusion of 
the test. 

Linear expansions of the three sandstones per- 
pendicular to the bedding are compared with the 
expansion of quartz in Fig. 2. Curves for expansion 
parallel to the bedding and cooling curves are not 
shown for Bandera and Boise because of their 
close similarity to the Berea curves. All sandstone 
samples showed permanent elongation upon cooling. 
The volume expansions shown in Fig. 3 were calcu- 
lated as the sum of the perpendicular linear expan- 
sion and two times the parallel expansion.° 


Interpretation of thermal expansion data is as 
follows. 


1. The quartz contents of the three sandstones 
ranged from about 50 per cent for Boise to 90 per 
cent for Berea, and yet the linear thermal expansions 
of the three sandstones were-nearly the same and 
about equal to that for quartz. The presence of 
quartz, probably above some minimum amount, 
appears to control the expansion behavior of sand- 
stones. 

2. Expansion of quartz in the direction perpen- 
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FIG. 2 — LINEAR EXPANSION OF SANDSTONES PER. 
PENDICULAR TO BEDDING. 


dicular to the C-axis has a predominating effect 
on the expansion characteristics of the sandstones, 
Since the amounts of expansion of the sandstones 
in the directions parallel and perpendicular to the 
bedding were approximately the same, random 
orientation of the quartz crystals would be expected. 

3. Permanent elongation of the samples after 
cooling resulted from deformation within the sample, 
due primarily to the differential expansion of the 
quartz grains. At temperatures above that for a-8 
quartz inversion, where the coefficient of expansion 
for quartz becomes negative, the role of other 
mineral constituents becomes important in determin- 
ing expansion and deformation characteristics of 
the sandstones. 


THERMAL REACTIONS 


In measuring thermal diffusivities of rocks by 
the unsteady-state method,? anomalies were ob- 
served to occur in the differential temperature 
records. Similar to differential thermal analyses 
(DTA), these temperature anomalies may be corre- 
lated with thermal reactions which occur in certain 
of the mineral constituents. The most consistent 
reaction is the inversion of quartz from the @ to 
B-phase at 573°C. The amount of heat needed to 
complete this inversion is known to be 4.825 cal- 
ories/gm.© The phase change is reversible, and an 
equivalent amount of heat is liberated upon cooling. 

Other typical thermal reactions which may occur 
are summarized in Table 1. Most of these reactions 
occur over a broad temperature range and require 
substantially larger amounts of heat than the quartz 
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FIG. 3 — VOLUME EXPANSION OF SANDSTONES. 
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TABLE 1 — HEATS OF REACTION FOR 
SEVERAL MINERALS?’ 


Heat of 
Temp. Range Reaction 
(°C) Mineral (Calories/gm) Reaction 





25-220 Ca-Montmorillonite 127 
25-220 Mg-Montmor illonite 135 


Desorption 
Desorption 


400-695 Mg-lllite 64 Decomposition 
455-642 Kaolinite 253 Decomposition 
554-723 Ca-Montmorillonite 67 Decomposition 

573 Quartz 4.82 a-B Inversion 
700-830 Calcium Carbonate 465 Decomposition 
790-950 Mg-lllite 15 Decomposition 


816-908 Ca-Montmorillonite 26 


Decomposition 





inversion. These reactions are normally irreversible 
or are only reversible under special conditions. 
Differential temperature records for the three sand- 
stones shown in Fig. 4 demonstrate this point. 
The dashed lines show the quartz inversion starting 
at approximately 575° C, followed by several other 
temperature anomalies. The solid lines show the 
thermal behavior of the same sandstone samples 
on a repeat run. The reversible quartz reaction is 
the only important anomaly remaining. 

The purpose of the present investigation was 
to estimate the magnitude of heat required to satisfy 
the thermal reactions and to compare this with the 
heat required to raise the temperature of the rock 
excluding thermal reactions. The known heat of 
reaction for quartz inversion and the known quartz 
content of Berea sandstone were used for standard- 
ization purposes. 

The apparatus used to evaluate thermal reactions 
is the same as that described earlier for thermal 
diffusivity measurements.? In the earlier work, 
cylindrical core samples were mounted in an electric 
furnace and heated at a constant rate of temperature 
rise. At temperatures where thermal reactions were 
known to occur, temperature was held constant 
until the reaction was completed. The edge temper- 
ature of the sample and edge-to-center differential 
temperatures were measured with thermocouples 
and traced on a recorder chart. It was shown 3 that, 
for these conditions, thermal diffusivity is inversely 
proportional to the differential temperature; i.e., 


ed 5 rr 


where D, is thermal diffusivity (D, = a in Refs. 1 
through 3), @ is the distance between center and 
edge thermocouples, 4 is the heating rate and 
AT is the temperature differential. 

In the present case, differential thermocouples 
were located at distances r/a = 1.0, 0.8, 0.6, 0.4 
and 0.2 from the center of the sample. The rate of 
temperature rise at the sample edge was maintained 
constant, and temperature differentials between the 
center of the sample and the several r/a locations 
were recorded. Fig. 5 shows the thermocouple 
responses caused by the a8 quartz inversion. 
Physically, the reaction is explained as follows. 
When the outer surface of the core reaches reaction 
temperature (573°C), part of the heat supply is 
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used to promote the reaction. The normal temperature 
gradient in the sample is disturbed and a slight 
decrease in the differential temperature is noted 
for the first thermocouple. As the reaction reaches 
and then passes the first thermocouple, a sharp 
increase in the temperature differential for this 
couple occurs. The edge temperature continues to 
increase at the constant input rate; but less heat 
reaches the unreacted portion of the sample because 
of the large amount of heat consumed by the reaction, 
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and the rise in center temperature lags. Thus, a 
sharp temperature discontinuity must exist at the 
reaction front. Similar response is noted as the 
reaction proceeds past each interior thermocouple. 
All curves show a maximum value when the reaction 
reaches the center thermocouple. The temperature 
differentials then decrease to re-establish the 
original temperature gradient in the sample. 

If the supply of heat to the system is constant, 
a heat balance may be written as 


(D.C, Sr) ~ (26, $2) =n $8 ee: 


where Cy is heat capacity, 57/dr is temperature 
gradient, Hp is heat of reaction, 5R/856 is the rate 
of movement of the reaction through the sample 
and a and f refer to the unreacted and reacted 
zones, respectively. Since the change of D, within 
the temperature range of the reaction is small, the 
assumption of constant heat supply is valid and 
Eq. 2 should apply to the present system. However, 
attempts to confirm the equation with experimental 
data from Fig. 5 failed. Many more thermocouples 
would be needed to establish reliable values of 
the temperature gradient ahead of and behind the 
reaction front, and the rate of movement of the 
reaction. 

Resort was made to an empirical method to obtain 
useful information from the experimental data. The 
method, which is used in quantitative interpretation 
of DTA data,’ involves the correlation of areas of 
the anomalies on the differential temperature-vs- 
edge temperature curves with heats of reaction. 
In Fig. 5 the areas between an estimated base line 
and the anomalies were measured for each thermo- 
couple location. These areas were then plotted 
against calculated heats of reaction for a volume 
of unit length of sandstone contained within each 
thermocouple radius and the center of the sample. 
A reasonably good straight-line relation resulted, 
as shown in Fig. 6. A value of 1 calorie/unit of 
anomaly was used in determining the heats of un- 
known reactions. 
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. 6 — AREA-HEAT OF REACTION CORRELATION. 


TABLE 2 — EXPERIMENTAL HEATS OF REACTION 


Temp. 
Range Heat of Reaction 
Sample (°C) 


Equiv. Heat 
(Calories/gm rock) (Calories/gm rock) 
577-615 2.9 VN 
615-735 21. 34 
735-750 -0.9 4 
750-825 5.6 22 
825-920 _ 79 27 

> 36.5 x 9 
565-600 4.3 9 
600-725 10. 34 
725-840 14, 33 
840-925 8.7 25 

> 37.0 > 101 
565-620 2.6 14 
620-725 3.9 28 
725-975 _ 71 





Bandera 





In analyzing the data of Fig. 4, the only known 
reactions were the a-$ quartz inversions. Differ- 
ences in the original and reacted curves below 
573°C are due primarily to the reduction in thermal 
conductivities of the reacted samples. This is 
probably caused by the loss of free and combined 
water and structural damage of the cores by differ- 
ential thermal expansion during the initial heating 
run. Reactions above 573°C include, in addition to 
quartz inversion, dehydroxylation of clay minerals 
and decomposition of calcium carbonate. It was 
not considered practical to attempt to separate 
specific mineral reactions because the purpose of 
the work was only to evaluate the additional heat 
requirements due to thermal reactions. Breaks in 
the curves were used to establish temperature 
ranges of the important reaction zones as shown in 
Table 2. Areas between the original and reacted 
curves of Fig. 4 were measured, and heats of 
reaction were calculated from the correlation factor 
and bulk densities of the samples. 

In Table 2, the heats of reaction are shown and 
compared with the amount of heat required to raise 
the temperature of the rock without thermal re- 
actions. The reaction heat is 27 per cent of the 
total heat requirements for both Bandera and Berea 
sandstones. A similar analysis was not made for 
Boise sandstone because of the difficulty of ex- 
plaining the large thermal reaction above 725°C. 
This same behavior persisted in tests on several 
other samples of Boise. 


SUMMARY AND APPLICATIONS 


The results of these tests, made at atmospheric 
pressure on oven-dried sandstone samples, typify 
behavior at just one set of limiting conditions. 
Before such data can be applied directly to sub- 
surface calculations, the effects of a wider range 
of conditions relating liquid saturation and over- 
burden pressure needed to be investigated. 

Differential thermal expansion of the mineral 
constituents of rocks will undoubtedly lead to 
structural damage under subsurface heating con- 
ditions. The magnitudes of effects on mechanical 
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and flow properties are unknown, but these currently 
are being investigated in this laboratory. It is 
possible that cores taken from high-temperature 
reservoirs may be altered by cooling as well as 
by the release of overburden and pore pressure. 

Thermal reactions involving the release of ad- 
sorbed and combined water and the liberation of 
such gases as carbon dioxide will be different 
under subsurface conditions of liquid saturation 
and pressure than the same reactions observed on 
oven-dried cores at atmospheric pressure. Tem- 
peratures of reaction will be higher and the heats 
of reaction may be different. Brindley and Nakahira® 
report that the starting temperature of the endother- 
mic kaolinite reaction is raised by 100°C under a 
water-vapor pressure of 6 atm. Further laboratory 
investigation will be necessary to evaluate these 
effects. 


NOMENCLATURE 
D, = thermal diffusivity (= ain Refs. 1 through 3) 


a = distance between center and edge thermo- 


couples 
b = heating rate 
AT = temperature differential 
r = radius of core sample 
Cy = heat Capacity 
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= temperature gradient 


oT 
or 
Hp = 
OR 
a0 

0 


heat of reaction 


rate of reaction movement through the sample 


time 


RD 
" 


unreacted and reacted zones, respectively 
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Three-Phase Imbibition Relative Permeability 


J. NAAR* 
JUNIOR MEMBER AIME 
R. J. WYGAL 


ABSTRACT 

An equation for three-phase (water, oil, gas) 
imbibition oil permeability is developed, assuming 
the water to be the dominant wetting fluid. Oil 
isoperms are obtained for consolidated sandstones 
characterized by 1/P,2 = CS*. The evolution of an 
oil-gas system imbibing water fromS,,; to(S., imb max 
is shown to proceed along a line of constant oil 
saturation with increasing oil permeability and 
decreasing gas saturations. When the gas saturation 
cannot be reduced further, the system evolves along 
a line of constant S, with decreasing oil saturation 
and permeability. The initial gas saturation is shown 
to reduce markedly the effect of complete wetting 
by either oil or water on flow performance. 


INTRODUCTION 


Imbibition oil isoperms are required for per- 
formance prediction when a well is producing water, 
oil and gas. This situation occurs in multiphase 
displacements such as underground combustion, 
steam injection and the water flooding of highly 
depleted reservoirs. 

In a recent paper,! a model was presented for the 
prediction of two-phase imbibition characteristics. 
This paper extends the imbibition model to the case 
of three phases by assuming that the water is the 
dominant wetting fluid. 

The following results were obtained from the 
model: (1) an analytical expression of oil isoperms; 
(2) oil isoperms as functions of reduced water, oil and 
gas saturations, valid for all sandstones having a 
capillary pressure curve which can be approximated 
by 1/P,* = CS*; and (3) evaluation of the three- 
phase flow performance as dictated by complete 
wetting by either oil or water. 

The agreement between predicted and observed 
oil recovery in the presence of a gas phase, reported 
in Ref, 1, is a partial support for the present develop- 
ment. However, experimental data are not available 
at this time to check fully the model predictions. 
Perhaps this paper will stimulate the collection of 
such data. 
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THEORETICAL 


The imbibition model of a porous medium has 
been described previously, and the reader is referred 
to the paper of Naar and Henderson! for details, 
In brief, the model is formed by the random inter- 
connection of straight capillaries, with a provision 
for the blocking of the non-wetting phase by the 
invading wetting fluid. 

The following expressions were derived from the 
model. 

1. The amount of oil blocked by the invading water 
is 


\-@ 
Sop = (1-5, ,) J xax - 
Oo 


where a is defined by 


1-Sy 


C=. * 
\- S,/ 


2. The imbibition relative permeability to the 
wetting phase is 


rnd /mb Ss; imo” 3" 
p. 2 


#3_,2 
k : -2 Ye ra 
riw,imb kA w,lmbdy 


. (18) 


The permeability k is written 


' e . 
k= gry? +2 as 
of 

3. Imbibition relative permeability to the nonwet- 

ting phase — for the nonwetting fluid during an 

imbibition process, the relative permeability at a 

given wetting-phase saturation S%, jmp is equal to 

the drainage relative permeability at a saturation 
S*y,drain defined by the relationship 


ia a 


* wT * | *#2 
Sw, imb ~ Sw, drain™ 2 Sw, drain - + (19) 
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THREE-PHASE IMBIBITION 
RELATIVE PERMEABILITY 


Relative Permeability to Water 


If oil and gas are considered equivalent to a single 
nonwetting phase, the system can be treated as 
having only two phases; the relative permeability 
to water is then equal to the wetting-phase perme- 
ability during drainage. If only the capillary pres- 
sure curve is known, the relative permeability to 
water can be computed from Eq. 18. 


Relative Permeability to Gas 


To compute the relative permeability to gas, the 
system may again be treated as a two-phase system 
with water and oil forming the single wetting phase. 
The permeability to gas is then computed as explained 
in Ref. 1. 


Relative Permeability to Oil 


The permeability to oil is more complex to com- 
pute. The mechanics of the displacement are de- 
scribed as follows and the analytical details are 
given in the Appendix. 

When the water invades the porous medium and 
occupies the pores of radius r which were initially 
full of oil, part of the oil in this size class of pores 
is blocked and the rest invades pores full of gas. 
In turn, part of the gas phase is blocked and the 
remainder is pushed out. The effects of these move- 
ments are twofold. 

1. Because a certain volume of oil is trapped, the 
oil permeability tends to decrease. 

2. Because a certain volume of oil moves from 
small pores into larger ones, the permeability tends 
to increase. 

The relative permeability to oil, therefore, is a 
balance between these two tendencies, and it is 
possible to observe an increase in oil permeability 
at constant oil saturation. 

After the oil invades all the size classes of gas- 
filled pores, the oil saturation and the oil relative 
permeability decrease when the water saturation 
increases. 

The permeability to oil at any stage of the imbi- 
bition process is computed by assuming that the 
oil phase is wetting with respect to the gas. The 
water saturation plus the blocked oil saturation 
form a fictitious, irreducible wetting-phase satura- 
tion. 

For consolidated sandstones having a capillary 
pressure curve approximated by the equation 1/P? = 
CS* kyo, imb iS given by 


“3% a 
Aro, imb = Sot (Sop + 3S,,) - °° °° (16) 
where S$, and S*;,, are defined by 


Ae ae 
aed Jf 
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Oil isoperms for consolidated sandstones are shown 
on a trilinear plot in Fig. 1 as functions of the 
reduced water, oil and gas saturations. The dis- 
placement mechanism proposed indicates that at 
the beginning of the imbibition process the water 
saturation S*¥, will increase at the expense of the 
gas saturation at constant oil saturation until no 
more gas is trapped. At this point the water satura- 
tion will increase at the expense of the oil satura- 
tion at constant gas saturation. In a trilinear dia- 
gram, therefore, the path describing the evolution 
of an oil-gas system is formed by two straight lines. 
By following the dotted lines from a point on the 
S*,, = 0 axis to the locus of the intersection of 
these lines, the history of any element during imbibi- 
tion may be traced. 


LIMITS OF WATER-OIL BEHAVIOR 
DICTATED BY WETTABILITY 


The influence of wetting on oil-water flow behavior 
has been examined in Ref. 1. The geometry of the 
porous medium, represented by S,,;, was shown to 
be important in assessing this influence. When gas 
is present as the third phase, the model and the 
assumed imbibition mechanism teach the following. 

1. The total oil recovery* will exceed the 50 per 
cent limit of the two-phase system and is given by 


J Sy/ 
Nef ips he eel: 

a ae Oe 

2. For equal values of the oil recovery, the 
value of k,/k,, is not a function of S,,; but is a 
function of the wettability and of the initial gas 
saturation. This dependence is shown in Fig. 3 
where (k,/k,,) water-wet/(k,/k,,) oil-wet is plotted 





* See Eq. 7 of Ref. 1. 





FIG. 1 — THREE-PHASE IMBIBITION. 
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FIG. 2 — THREE-PHASE DRAINAGE. 


vs Sg; for N,/N = .40. This ratio decreases rapidly 
when S,; increases. From the standpoint of oil re- 
covery and flowing water-oil ratio, the higher the 
initial gas saturation, the less significant the wet- 
tability of the medium becomes. The values of 
ko/k for an oil-wet medium and for increasing 
water saturation were computed from Eqs. 23 and 
24 of Ref. 1. The water isoperms for that case are 
shown in Fig. 2. 

3. The water saturation at a given recovery is a 
function of S,,; and the initial gas saturation. This 
can be seen in Fig. 3 where (S,, ) water-wet/(S,,) oil- 
wet is plotted against S,; for S,,; = .15 and S,,; = 
30. This ratio increases with S,; and its rate of 
increase is an increasing function of S,,;- The 
physical significance of this increase is immediately 
evident — water must be imbibed in sufficient 
quantity to displace the gas phase before the oil 
phase can be expelled. 

The proposed displacement mechanism applies 
to a small elemental volume (AV) of porous rock. 
The gross production behavior then can have the 
appearance of an initial gas production followed by 
an oil bank. This type of apparent behavior was 
observed by Kyte, et al.3 

The observation that, during imbibition, the gas 
phase is produced first was also made in this 
laboratory. We are indebted to J. J. Taber for the 
curves of Fig. 4 which show the imbibition history 
of a Berea core. The agreement between observed 
and predicted final gas and oil recoveries reported 
in Ref. 1 further supports this theory. It is also 
felt that, while experimental oil isoperms are not 
available at this time, the success of the predictions 
of relative permeability and oil recovery in two-phase 
flow are strong points in favor of the present develop- 
ment. 


CONCLUSIONS 


A mechanism to describe three-phase water 
imbibition has been proposed and applied to a 
theoretical model of a consolidated porous medium. 
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FIG, 3 — INFLUENCE OF WETTABILITY ON FLOWING 
OIL-WATER RATIO AND WATER SATURATION AT 40 
PER CENT RECOVERY. 


A formula has been derived for the imbibition oil 
permeability and has been used to compute oil iso- 
perms for consolidated sandstones. 

The gas phase has been shown to decrease the 
influence of wettability and to be important in 
determining the amount of water required to effect 
a given oil recovery. 


» NOMENCLATURE 


constant 

capillary pressure 
saturation 

total pore volume 
pore volume produced 


ratio of the pore volume full of nonwetting 
fluid to the total pore volume 


permeability 

pore radius 

porosity 

interfacial tension 
SUBSCRIPTS 


drain = drainage 


imb = imbibition 
g= gas 
o = oil 
oil blocked 
= oil flowing 


relative 
water 
= water irreducible 
fictitious water irreducible 
SUPERSCRIPTS 
* indicates a reduced value of saturation or 


; 
i.e., N Ses = wv we (see 


porosity; > 


Appendix) 
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FIG. 4 — IMBIBITION OF WATER INTO A BEREA CORE CONTAINING WATER, OIL AND GAS. 
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APPENDIX 
The fictitious, irreducible wetting saturation 
which is needed to compute oil isoperms requires 
the definition of new symbols to avoid confusion. 
The final results, however, are expressed in terms 
of the reduced saturations S%,, S*%, and S*,. The 


new symbols are the following. 
1. Oil blocked by the invading water, 


Sw Swi 


1- Sy 
Sop = (I-S,; WA LA ate o &@ * (1) 


2. Oil saturation free to move, 
= = oe ee ee ee 
Sor = Sy - Sop (2) 
3, Fictitious, irreducible water saturation, 
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Sty = 5, + Sy, a: ee owecletee s See 


4. Reduced saturations (with respect to S,,;), 


oe 9 Say 9s 4 eS Ce 
| -S,, 
eee Ss 
1-S,, 
* So 
—. Vere. 
0 1- Swi 
Spe ~ Sui 


2 Piet ow eoeeesanle 7 
i xy m7 


5. Reduced saturations (with respect to Sy,,), 


Pt 
| - Sy, 
S 
Pt of 
= ee . . . 9 
Se * ae (9) 


From Eqs. 7 and 8 and from Eqs. 7 and 9, the 
following relationships are obtained between S* and 
Se, Sor and S*3s . 












S*- Sh, 


als = eee laa . . . o . . . . . (10) 
i- Sz, 
Sor 

-_ Se lw ee é& -& s .© 7 7 . e (11) 

al \- Sty 
In addition, S$, is given by 
“2 
a ae (12) 


COMPUTATION OF THE OIL PERMEABILITY 


Consider the medium at irreducible water satura- 
tion S,,;, oil saturation S$, and gas saturation S 
and let the water invade the medium to reach the 
value S,,. A part of the oil is blocked by the water, 
(S,,); however, the gas phase is unable to dis- 
tinguish between the water phase and the blocked 
oil. Only the oil free to flow contributes to the 
oil permeability; for all practical purposes, the 
water and the trapped oil can be considered as 
immobile and equivalent to a fictitious, irreducible 
wetting-phase saturation. The permeability to oil, 
therefore, will be that of a medium having an irre- 
ducible wetting-phase saturation equal to S,, + So, 
and a wetting-phase saturation free to flow equal 
to S, — Sop. It is equal to 


* 3 
bo, imo =| 8 -Syy)] Y x Sor 


I oe a ——" 
ti 4 5 as. 
(7) °, ‘ 





- (13) 


(See Ref. 1, Eq. 30.) The relative permeability to 


oil is obtained by dividing ko jm, by &, which is 


given by 








t * 
KP U- Sy) 7? f ES gst. (14) 
o FP 
Then kyo imp is written, 
3 
_ !- Sty Ht 
kro, imb * er Sof 
Sy x -s" il 
ya oe 
(Z c (15) 
| 1-s* ie 
f es 
oO P, 


Oil isoperms can be calculated from Eq. 15 for any 
consolidated porous rock. 

For consolidated sandstones, the capillary pres- 
sure curve can be fitted by 


* 


ee ee ee ee ee 


: 7 

p,* 

Replacing in Eq. 15 and integrating gives 
* 


3 * * 
hin sme * Gy th, + 88): > 


At the beginning of an imbibition process, Sj, = 0, 
S%b = Oand kyo, imb = 554.. This value is in agree- 
ment with the results of Wyllie and Gardner (Eq. 
31 of Ref. 2). * kk 
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Correlation of Surface and Interfacial Tension of Light 


Hydrocarbons in the Critical Region 
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ABSTRACT 


Empirical equations for surface tension of propane 
and normal butane as functions of reduced temper- 
ature are obtained from experimental data. Another 
correlation relating surface tension to enthalpy of 
vaporization is given for these two compounds. In 
addition, new parachor numbers are calculated for 
the normal paraffin hydrocarbons. These numbers 
are utilized for the calculation of interfacial tension 
of two-component systems as functions of pressure 
and temperature, using a modified form of Weinaug- 
Katz equation. The experimental data for two binary 
systems are approximated by the correlation. From 
these results it is found that the interfacial tension 
in the high-pressure region remains extremely low 
at large pressure decrements below the critical 
pressure. Thus, it appears that condensate systems 
may have flow characteristics almost like single- 
phase conditions even though the pressure is within 
the two-phase region. 

Experimental data have shown that interfacial 
tension divided by density difference approaches 
zero as the critical pressure is approached. A 
calculation of wetting-phase saturations indicates 
that the saturation gradient at the two-phase contact 
becomes increasingly abrupt as the critical pressure 
is approached. 


DISCUSSION 


Prediction of the surface and interfacial tension 
of the light hydrocarbons and of two-component 
hydrocarbon mixtures at various temperatures and 
pressures may be made if other physical properties 
are known. Extensive experimental work on single- 
component and binary systems! is the basis for 
the correlations outlined in this paper. 
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Interfacial tension is defined as the specific 
surface-free energy between two phases of unlike 
fractional composition, while surface tension is 
defined as the specific surface-free energy between 
two phases of the same fractional composition. The 
usual definitions relating interfacial tension to a 
liquid-liquid interface and surface tension to a 
gas-liquid interface are not clearly defined when 
the critical region is included, and there is no sharp 
distinction between a gas and a liquid phase. 

Interfacial tension is probably the most important 
single force that makes one-half to one-third of the 
total oil actually in place in a reservoir rock un- 
recoverable by conventional gas-drive or waterflood 
methods. A rough estimate of this figure for the 
United States is 100 billion bbl. Interfacial tension 
presently is used by petroleum engineers in the 
estimation of saturation gradients at the gas-oil 
contact and at the oil-water contact. The data in 
this paper should prove useful for estimates of 
reserves involving gas-oil contacts. Relative per- 
meability undoubtedly is influenced by interfacial 
tension, for sufficiently small values. These data 
should be useful in determining how small the 
values are. In addition, these data should eventually 
add to our fundamental knowledge of surfaces. At 
the critical point, all surface excesses approach 
zero and the thickness becomes very large. 


SINGLE-COMPONENT SYSTEMS 


It has been observed! that the following relation- 
ships are good approximations to the physical 
properties of propane and n-butane. 


For propane, Ap = 0.80 (1 es eee? y = 49.5 
(1-~T/T,)1-20, y/Ap = 61.9 (1—T/T,)9-874, and 
y = 112.5 Ap3-68, 

For n-butane, Ap = 0.86 (1—T/T,)®-333,y = 52.5 
(1-T/T,)1:22, y/Ap = 61.0 (1-T/T,.)%887 and 
y = 91.2 Ap3-66, 

Guggenheim’s? values for these constants, not 
specifically for hydrocarbons, are Ap = C’ (1- 
T/T,)!/3, y= CC” (1-T/T,)!1/9, y/Ap = C**" (1- 
T/T,.)8/9, and y = CAp11/3, 








The symbols and units are 
Ap = difference in density of the liquid 
and vapor phases, gm/cc, 


y = surface tension, dynes/cm, 
T = absolute temperature, °R or °K, 
Te = critical temperature (absolute), °R 


or °K, and 
C4C% C’”’ = constants. 

The exponents for propane and n-butane are taken 
from log-log plots. The last equation may be recog- 
nized as the surface tension-density difference re- 
lationship leading to Sugden’s parachor,> with the 
exponent ‘°3/11’’ rather than ‘‘1/4’’. That is, 

pewom, ae 


\p 


where P =parachor, and 
M =molecular weight. 

From this modified equation, new values of 
parachor were calculated for some of the normal 
paraffin hydrocarbons. Surface-tension and density- 
difference data for methane, ethane, n-decane and 
n-eicosane were taken from Rossini‘ for n-pentane 
from Sage and LaceyS and from Katz,6 and for pro- 
pane and n-butane from Sage and Lacey and from 
experimental work! The values are shown in Table 
1. Because the parachor is (ideally) independent 
of temperature, the values of surface tension and 
density difference were selected at the temperature 
at which the best data were obtainable or at the 
temperature at which the vapor density is negligible 
if this latter data were not available. Fig. 1 shows 
a regular increase in the parachor with an increase 
in molecular size from 3 to 20 carbon atoms. Methane 





TABLE 1— PARACHOR OF 
NORMAL PARAFFIN HYDROCARBONS 





Density Surf. or 
Carbon Temp. Diff.  Intf. Tens. 

Compound Atoms Cr) (gm/cc) (dynes/em) Parachor 
Methane 1 -256 0.4218 13,78 77.9 
Ethane 2 -238 0.620% 26.34% 118.0 
Propane 3 100 0.443 5 5.47 158.1 

5.45% 158 
n-Butane 4 120 0.5325 9.27 200.5 

9.20 200 
n-Pentane 5 100 0.6035 14.0% 246 
n-Decane 10 68 0.730R  23,92R 463 
n-Eicosane 20 100 0.7754 = 27,2R 899 


R=Rossini (Ref. 4); S=Sage and Lacey (Ref. 5); and K=Katz 
(Ref. 6). 





TABLE 2 — ENTHALPY OF VAPORIZATION AT 
CONSTANT PRESSURE — PROPANE 


Enthalpy of 











Temperature Vaporization* Surface Tension** 
(°F) (calories/gm-mol) (dynes/cm) 
100 3271 5.48 
130 2878 3.66 
160 2366 1.97 
190 1660 0.58 


*Sage and Lacey (Ref. 5). 
**Current paper. 
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and ethane appear to diverge slightly from this 
pattern. Some of this divergence is probably caused 
by the measurement of surface tension against air 


rather 


vapor density. ; 

A correlation suggested by Magaril’ relates the 
enthalpy of vaporization at constant pressure to the 
surface tension. Tables 2 and 3 and Figs. 2 and 3 
express this relation for propane and for n-butane, 
The straight-line relationship on log-log paper 
appears to deviate slightly as the critical point is 
approached. 


TWO-COMPONENT SYSTEMS 


Prediction of interfacial tension in two-component 
systems is much more difficult than in the single- 
component systems because of the added composition 
variable. The correlation suggested by Weinaug and 


Katz8 


slightly modified for use with molal volumes rather 
than specific gravity. The relation is 


us | 
y as [Py (x4/Vx — ¥1/Vy) 


where 


1000 


900 


800 


P, = parachor of the first component, 
Py = parachor of the second component, 


Ty yy TT 









than its own vapor, and by the neglect of the 












has been used with the new exponent and 


k 
+ P2(x9/Vy - a/v, )]) 


y = interfacial tension, dynes/cm, 
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TABLE 3 — ENTHALPY OF VAPORIZATION AT 
CONSTANT PRESSURE — NORMAL BUTANE 








Enthalpy of 
Temperature Vaporization* Surface Tension** 
(°F) (calories/gm-mol) (dynes/cm) 
100 4816 10.49 
130 4566 8.67 
160 4266 6.88 
190 3912 5.14 
220 3498 3.48 
250 2980 2.02 
280 2307 0.80 


*Sage and Lacey (Ref. 5). 


**Current paper. 





xy, X2 = mol fraction in the liquid phase, 


yj, ¥2 = mol fraction in the vapor phase, 
vy Vy = molal volume, cu ft/lb mol, and 
k = new exponent of the parachor relation 


= 11/3. 

The modified values of the parachor together 
with the new exponent are used in this equation 
to predict values of interfacial tension determined 
experimentally! on the methane-normal pentane 
system and on the methane-normal decane system. 
These experimental data are summarized on the 
pressure-temperature plots in Figs. 4 and 5. The 
boundaries of this region are the vapor-pressure 
curves for the pure components and the line of 
critical loci of the various intermediate mixtures. 
Superposed on these diagrams are the lines of con- 
stant interfacial tension. The known data for these 
iso-interfacial tension lines are the surface-tension 
data of the pure compounds and the experimental 
data by the authors, which include a very limited 
range in comparison to the complete two-phase 
region. Since the properties of the liquid and the 
gas at the critical point are identical, the: line of 
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critical loci is also the line of zero interfacial 
tension; hence, the subsequent iso-interfacial 
tension lines must roughly parallel this boundary. 
These data and assumptions allow extrapolation 
and interpolation to all points within the two-phase 
region. It is evident from these plots that the inter- 
facial tension of a binary mixture in the temperature 
range in which the critical pressure is high increases 
much more slowly with a decrease in pressure than 
either of the pure components. In the case of the 
methane-normal pentane system, a pressure drop 
of 450 psi from the critical pressure is required 
at 100°F before the interfacial tension reaches 1 
dyne/cm. The same increase for the pure components 
occurs with a pressure drop of about 120 psi. In 
the methane-normal decane system the effect is 
even more pronounced, the pressure drop at 100°F 
being 1,040 psi for a comparable increase in inter- 
facial tension. These facts suggest that interfacial 
forces may be negligible, as far as flow conditions 
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are concemed, for a condensate reservoir for an 
appreciable portion of the two-phase region. 

Prediction of these experimental results from the 
modified equation of Weinaug and Katz was attempted 
with the PVT data of Sage and Lacey. A summary 
of these calculations is given in Table 4. 

The equation for binary systems may be generalized 
to multicomponents systems as follows. 


oa 9 (F =) P{ 22 
y* 62.43 [* W, ~v,/* Av, ~V 


y 
p (22 m\{( 
+ P,, Vz Vy . 


This calculation requires a knowledge of the 
composition of both phases and the molal volume 
or density of each phase. In the absence of density 
data, reduced temperature calculations may be used 
as an approximation to the density. The parachor 
numbers for each compound of the normal paraffin 
hydrocarbons from methane to n-eicosane are shown 
in Fig. 1. 


SURFACE EXCESS 


Excess concentration of methane at the liquid 
surface for the two-component systems may be ap- 
proximated from the equation given by Rice.9 


I; = -pn (x) 
T 


where I’; = excess of methane in the surface, mol/ 
sq cm, 


Pn = density of methane, mol/cc, 
p = pressure, dynes/sq cm, and 
T = temperature, °R or °K. 


It must be noted that this equation was derived 
with the assumption of a zero surface excess of the 
other (heavier) component. 

Calculation of surface excess from experimental 





TABLE 4 — INTERFACIAL TENSION CORRELATION 
(Weinaug and Katz®) 


Interfacial Tension 














Temperature Pressure (dynes/cm) 

Components rr (psia) Exper. Cale. 
Methane and 100 600 9.02 9.30 
n-Pentane 1250 4.59 4.78 
1750 2.01 2.09 
2250 0.27 0.21 
220 600 4.62* 5.02 
1250 1.98* 2.14 
2000 0.07* 0.05 
340 800 0.37* 0.37 

Methane and 100 1250 11.3 11.8 
n-Decane 2000 7.35 7.02 
3500 2.40 1.84 
5000 0.16 0.09 
280 1000 9.13* 9.30 
3000 2.31* 2.32 
4500 0.06* 0.02 
460 1000 3.3* 3.95 
2750 0.3* 0.10 

*Extrapolated 
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values is obtained by using the slope of isotherms of 
interfacial tension plotted against pressure for both 
two-component systems. If pressure is expressed 
in pounds per square inch, interfacial tension in 
dynes per centimeter and specific volume in cubic 
feet per pound-mole, a conversion factor must be 
used in the modified equation. 


r, — — 14.013 x 10! / dy 
ai Vi Op }, 


For example, methane-normal pentane at 100°F 
and 1,500 psia, 





Ce) 

52) = -0.00506 a. 

P/ 
V; = 3.50 cu ft/lb mol (methane), and 
I, = 2.03 x 10!4 mol methane/sq cm. 


Results of surface-excess calculations are given 
in Tables 5 and 6. As the pressure approaches 
zero, the value of the density term in the surface- 
excess equation approaches zero and the isothermal 
change in interfacial tension with pressure remains 
finite; hence, the surface excess as the pressure 
approaches zero also approaches zero. See Figs. 6 
and 7. As the pressure approaches the critical point 
of the mixture, the density of methane remains finite; 
but, according to experimental results, the change 
in interfacial tension with pressure approaches 
zero. Therefore, the surface excess again approaches 
zero as the critical pressure is approached. For 
intermediate values of pressure, the density term 
is positive and the slope was found to be negative; 
hence, the surface excess is positive between zero 
and the critical pressure, and has at least one 





TABLE 5 — SURFACE EXCESS OF METHANE IN THE 
METHANE-NORMAL PENTANE SYSTEM 
(MOLECULES x 10°!4/sQ cm) 


Temperature (° F) 





Pressure 
(esi) «=» «0 150.———NKD.7 200 
2250 1.38 0.72 0.57 0.42 - 
2000 1.90 1.34 1.25 1.17 - 
1750 2.09 1.53 1.40 1.36 0.99 
1500 2.03 1.50 1.40 1.33 1.03 
1250 1.92 1.46 1.35 1.27 0.98 
1000 1.66 1.26 1.16 1.08 0.90 
800 1.40 1.07 0.98 0.93 0.76 
600 1.17 0.82 0.78 0.74 0.60 





TABLE 6 — SURFACE EXCESS OF METHANE, 
METHANE-NORMAL DECANE (MOLECULES x 10°!4/sQ CM) 


Temperature (°F) 





Pressure 
(psia) 100. 130 160 
5000 0.98 0.86 0.69 
4500 1.42 1.39 1.33 
4000 1.74 1.56 1.57 
3500 2.03 1.93 1.83 
3000 2.34 2.21 2.12 
2750 2.45 2.38 2.20 
2500 2.54 2.41 2.22 
2250 2.54 2.30 2.13 
2000 2.44 2.12 1.96 
1750 2.29 2.04 1.88 
1500 2.00 1.79 1.67 
1250 1.66 1.50 1.40 
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maximum point at some intermediate pressure. 

The occurrence of a maximum, followed by decreas- 
ing values of surface excess, is attributed to the 
increased value of the surface excess of the heavier 
component. This value is probably small in the 
low-pressure region, but it becomes large as the 
pressure becomes large. Thus, it is reasonable to 
assume that the decrease in the ‘surface excess’’, 
as calculated, is the result of the surface excess 
of the heavier component’s becoming significantly 
large. The maximum value of the calculated surface 
excess of the methane-normal pentane system (or the 
methane-normal decane system) might be associated 
with the attainment of a monolayer of methane on a 
solid surface, since extrapolation of these values 
for the various isotherms gives fairly good agree- 
ment of values of surface density!9 of a monolayer 
of methane on a solid surface at ~297°F. 

It appears that surface excess might be useful 
as a correlating device for interfacial tension, since 
its values appear to be much more constant for the 
various systems than interfacial tension itself. For 
instance, a plot of surface excess divided by density 
difference appears to give nearly a straight line 
when plotted against pressure and might possibly 
yield a finite value at the critical pressure, although 
in the case of the methane-normal pentane system 
the value appeared to approach zero very near the 
critical point. 

Experimental datal have shown that not only the 
interfacial tension, but also the interfacial tension 


T 
| 
= 2.5210"|— | —— $9} + 





} 

2.0x10* 
| 
(oe 
| 











SURFACE EXCESS [ motecuies « cm’ 





A I A 
O 500 1000 1500 2000 2500 3000 3500 4000 
PRESSURE [opsie) 


FIG. 6 — SURFACE EXCESS, METHANE-NORMAL PEN- 
TANE. 


25,10" 
20.10: 
15210") 
10110"; 


05:10") 





SURFACE EXCESS [ motecutes x em 


o 


© 1000 2000 3000 4000 5000 6000 7000 8000 
PRESSURE [psie] 


FIG. 7 — SURFACE EXCESS, METHANE-NORMAL 
DECANE. 


DECEMBER, 1961 





divided by density difference approaches zero as 
the critical point is approached. The significance 
of the latter is emphasized by a calculation of 
wetting-phase saturation of porous media from inter- 
facial-tension data. Plotting saturation gradient 
vs height, one obtains a very sharp gas-liquid 
contact as the critical point is approached! This 
sharp liquid-gas saturation gradient may be a con- 
tributing factor to the nonequilibrium conditions 
found in oil reservoirs by Cupps, Lipstate and Fry.}! 


SUMMARY 


A considerable amount of data has been accumu- 
lated on the surface tension of propane and normal 
butane, especially near their critical temperatures. 
In addition, the binary systems — methane-normal 
pentane and methane-normal decane — have been 
investigated throughout a wide pressure and temper- 
ature range, particularly in the region near the 
critical points. All these data point to the regularity 
of the pressure and temperature functions of surface 
and interfacial tension. Of most interest in the 
data, however, is the persistence of low interfacial 
tensions with relatively large changes in pressure 
in the region about the critical points. The importance 
of this fact to the petroleum engineer lies in its 
relationship to liquid saturations and flow behavior 
in an oil reservoir. 

A means of predicting interfacial tension for 
hydrocarbon mixtures from composition and volumetric 
data appears to fit experimental results very well. 
Other correlating devices such as reduced temper- 
ature relations, enthalpy of vaporization relations, 
and surface-excess calculations give promise of 
relating the surface phenomena to better known 
physical phenomena. 
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An Approximate Method for Computing Nonsteady- 
State Flow of Gases in Porous Media 


L. G. JONES 


JUNIOR MEMBER AIME .... 


ABSTRACT 


An approximate method of calculation is developed 
in this paper which allows duplication of radial 
unsteady-state gas flow computer results where 
Darcy’s law applies, such as those reported by 
Aronofsky and Jenkins and Bruce, Peaceman, 
Rachford, and Rice. Moreover, the new calculation 
method can be used to obtain results for radial 
unsteady-state gas flow obeying the quadratic flow 
law proposed by Duwez and Green. 

Means are discussed for predicting well behavior 
at single or superimposed flow rates in finite or 
infinite reservoirs, determining reservoir rock prop- 
erties from well-test data, reproducing and inter- 
preting back-pressure test data, and determining 
the radial extent and reserves of gas reservoirs 
from well-test data. Example calculations are 
presented for gas flow following both Darcy’s law 
and the quadratic flow law. 


INTRODUCTION 


Since the publishing of U. S. Bureau of Mines 
Monograph 7,1 most gas-well testing methods have 
been based on the equation g = y (p77 ~ p¥*, 
where q = production rate, py and p,, are formation 
pressure and sandface pressure, respectively, and 
y and a are constants to be obtained from test data. 
These methods, used for predicting both deliverability 
and ‘‘open-flow’’ capacity of gas wells, have been 
useful and accurate in many cases but unsatisfactory 
in others. Even at best, however, they do not supply 
information about the formation or lead to an under- 
standing of nonsteady-state gas flow in porous 
media. 

Many theoretically based studies2-4 of gas flow 
obeying Darcy’s law have been made.* Since the 
partial differential equations which result from 
combining Darcy’s flow law with the continuity 
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*These are examples, there are many others. 


264 


SOCONY MOBIL OIL CO., INC. 
DALLAS, TEX, 


equation are nonlinear, most of the published research 
consists of either numerical solutions or analytical 
solutions for linear approximating equations. Such 
solutions have been of limited value in field work 
due to their unhandy form and their failure to correlate 
most field data. 

There is evidence5,6 which indicates that Darcy’s 
law is inadequate to describe gas flow at some flow 
rates of practical interest. A quadratic flow law, 
which reduces to Darcy’s law at low rates, is more 
successful in accounting for experimentally observed 
behavior. This flow law has been applied success- 
fully to a few hypothetical reservoir cases in work 
which has not yet been published. However, these 
numerical solutions of the equations involved have 
been successful only on a special analog computer. 
Routine applications to field cases would be awkward 
and have not been attempted. 

The present paper describes an approximate 
method for computing nonsteady-state gas flow 
solutions which has been completely successful 
in predicting the results for both Darcy flow and 
quadratic flow obtained by elaborate numerical 
methods. The new calculation method allows deter- 
mination of the observable variables in gas-well 
testing at constant rates. It is similar to the scheme 
of using a succession of steady states suggested 
by Muskat7 in that it makes use of steady-state 
and material-balance equations. It also is similar 
to the work of Aronofsky and Jenkins? in that the 
new method includes Aronofsky and Jenkins’ stabil- 
ized flow equation as a special case. It improves 
upon both of these calculation schemes in that it 
accurately describes all portions of reservoir history 
and suggests means of determining reservoir rock 
properties from well-test data. 


This paper deals only with production from a 
reservoir, in which case the rate is defined as being 
negative. The reservoir studied here is a homoge- 
neous, disk-shaped porous body of uniform thickness, 
with all boundaries sealed except the inner radial 
boundary of the well. 
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THE CALCULATION METHOD 
ASSUMPTIONS AND DEFINITIONS 


The method is based upon the following assump- 
tions and definitions. 

1. The drainage radius rg is defined from integra- 
tion of the radial flow equation until this equation 
indicates that ry =r,, the radius of the outer extremity 
drained by a well, after which rg is defined as being 
equal to fo. 

2. All production is assumed to have come from 
within the drainage area bounded by rg as defined 
herein, and all of the reservoir energy is assumed 
to be provided by expansion of the reservoir gas. 

3. It is assumed that the average pressure Paye 
within the drainage area occurs at a characteristic 
radius fgye, which can be described by the relation 
fave/td = C, where C can be uniquely determined 
as a function of 1y/re. 

4, The “‘early transient period’’ is defined as 
that period during which rg < re. The “‘post early 
transient period’’ is defined as that period during 
which rg = fj. 


DARCY FLOW 


The radial flow formula for gas flow obeying 
Darcy’s law may be written as 


7 kbMp 2 
os sie tampa ee me a 
In rp uWzaeRT (Pp Pwo*), 


or 


In rp = (p p2 — Pwo), a @. 6) eee (OR 


WD Zave 


—-W uRT 


where Wp = dimensionless mass flowrate, ___' -_, 


Tp = dimensionless radius, r/r,,, 


Pp = dimensionless pressure, p/p;, 

k = permeability, 

hb = net pay thickness, 

M = molecular weight, 

W = mass flow rate, pq, 

ft = viscosity, 

Zave = average compressibility factor, 
R = gas constant, 
T = absolute temperature, 

bwp = dimensionless well pressure, 


pb; = initial pressure, 


q = volumetric flow rate which is defined 
as being negative for production, 


p = density, and 
Ty = well radius. 
Early Transient Period 


During the early transient period the pressure at 
the drainage radius is defined as being equal to the 
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shut-in or initial pressure p;. The drainage radius 
for this case is then defined by the equation 


1 
Inrgp = Wien (pip? - Pwo?) . «++ (2) 


where rgp = dimensionless drainage radius and 
Pp = dimensionless initial pressure. 


Since 
C = Taye/Td: 
Td =TdD'w 
and 


Tave = TaveD'w 


C = faved /TaD 


or 
fue @ Chap. +s 6 ae eee elec e og 
By definition, Pgye occurs at rgye and 
In TayeD = Wpzave (PaveD” — Pwo?) 
= InCrap = Inc + Inrgp » a ae (4) 


where P,,,op is dimensionless average pressure. 


The average pressure decline, Apgye, within rgp 
is defined as p; — Pgye- This pressure decline must 
account for total production G,,, in moles, so that 


ZquaGekt 


»* » 
bdmrg? -r 2) ” 


\P ave 


where ¢ is the porosity of the reservoir. One may 
note the relationship between production and rg. If 
APave is constant, G,, is approximately proportional 
to rqg%. This indicates that rg increases very rapidly 
upon initiation of flow, and very slowly after rg 
becomes. large. 

Another bit of information can be gained at this 
point by adding the negative of Eq. 2 to Eq. 4. The 
sum of these is 


T aveD se. 
042) - tn ”~ a. (1~Pavep2) . (6) 
It may be seen from this equation that, if C is 
constant for a constant-rate case during the early 
transient period, Pgyep and Paye must also be 
constant. 

Numerical digital-computer solutions for unsteady- 
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state flow, similar to those of Aronofsky and Jenkins2 
or Bruce, Peaceman, Rachford and Rice,3 were used 
for determination of the C values. In this case the 
results were obtained by using a central-differencing 
scheme to numerically solve the equation. 








0p p? = e2% op 
@u2 ap 






D;t 


2° 
Tw 








where u = In rp and tp = Only boundary 






conditions, values of p,, and production G,, must 
be known to détermine C. 

C values for the early transient period can be 
calculated by using Eqs. 2,4 and 5 in the following 
manner. 

Step 1 — Calculate rgp using Eq. 2. 

Step 2 — Using rgp from Step 1, calculate Apaye 
using Eq. 5. Then paye =P; — APave- 

Step 3 — Using Paye from Step 2, calculate C 
using Eq. 4. 

Fig. 1 shows values of C calculated from various 
computer runs. Rate is given in dimensionless 
terms for reasons explained in the Appendix. It is 
apparent that C is a function of yp /tep only. Two 
points are appreciably distant from the curve. One 
of these is for a dimensionless rate of .001, which 
is so low that digital round-off errors are appreciable. 
The other is for an rgp value of around 10. Since C 
= TayeD/Tdp it is apparent that, if C < 1, rayep 
will be < 1 when rgp = 1. Since this is physically 
impossible, the correlation obviously cannot work 
for extremely small values of rgp. This situation 
is not of practical importance, however, since it 
would prevail for only a few seconds in a well test. 











































Post Early Transient Period 


During this period rg = rg, and pe, the pressure 
at the outer boundary, declines. For this situation, 


In tr agvep = In Crep = InC +In rep 
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ZaveGumKT Z ave G»RT 
hd nm (re2 — 1,2) bd a re? 
ee ee + . . . * e es . . > = Ss (8) 





Abave = 


C values for the post early transient period can be 
calculated by using Eq. 8 to calculate pgye and by 
using this in Eq. 7 to determine C. Some calculated 
values for C in the post early transient period are 
presented in Table 1. In the table, Gp is dimension- 
less cumulative production and is the ratio of total 
production to total moles initially present in the 
reservoir. Values vary within the range .4486 to 
4758. The variation is probably due to a relaxing 
of the convergence limits in the digital runs. It is 
believed that a C value of .470 will adequately 
represent the entire post early transient period. This 
compares with a value of .472 found by Aronofsky 
and Jenkins. The difference is of no consequence 
and is due to the normal variations in digital soly- 
tions. 


Calculation Procedures 


If C, the gas properties and the reservoir con- 
stants are known, reservoir performance at constant 
flow rates can be predicted by using Eqs. 2, 3, 4 
and 5 for early transient periods, and Eqs. 7 and 8 
for post early transient periods. 

For the early transient period, the calculation 
procedure is as follows. 

Step 1 — Assume a value of p,,. 

Step 2 — Calculate ry using Eq. 2. 

Step 3 — Calculate ryp/rep; read C from Fig. 1. 

Step 4 — Using C from Step 3, calculate raye 
using Eq. 3. 

Step 5 — Using raye from Step 4, calculate paye 
using Eq. 4. Then Afgye = 2; — Pave: 

Step 6 — Using Apaye from Step 5 and rg from 
Step 2, calculate G,, using Eq. 5. 

For the post early transient period, the procedure 
is the following. 

Step 1 — Assume a value of py, 

Step 2 — Using C = .470 and p,, from Step 1, 
calculate paye using Eq. 7. Then Apaye = Pj — Pave 

Step 3 — Using Apgye from Step 2, calculate G» 
using Eq. 8. 











TABLE 1 
Ge c Gp c 
Run 103 Run 105 
.008568 4692 .008994 4517 
.009486 4613 .01470 .4695 
.009742 4657 01964 4486 
.01257 .4573 .03394 .4734 

01470 4638 

01753 4569 Run 106 
.01966 .4633 009459 .4723 
.03240 4595 .01802 .4758 
.05832 4579 .03037 4698 
.08369 4635 .05991 4659 
12285 4611 

14298 475) 
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QUADRATIC FLOW 


To adequately describe gas flow at moderate or 
high rates, it is necessary to add another term to 
Darcy’s flow law, 





where B is a parameter to be determined experimen- 
tally by a laboratory procedure somewhat similar to 
that for determining permeability. Multiplying this 
by p gives 





-p dp _ ue q Bp? q? 
dr k A AZ 
which equals 
Ww w2 
= B ere st 





+ 
k 2nrh 4 m2r2p2 


. . M . - 
Substituting rp =1r/r,, and p a Eq. 10 
ZeugRT 
and rearranging gives 











Ma dp? yw Ww Bw2 
iy Na eee ap = .(11) 
ZevehT drp k Tpb 27rp Tw? 
If 
7 ae 
a 
_Wp 
ins kh 
and 
w2 
)-—£,, Mwave ees: ER 
2 71h 
then 
pwr tw ee ee ee ee OD 


drp Tp Tp 


Rearranging variables and integrating yields 
¥ J fz 
p2 . p2,, a lnrp -7 (* _- ) 4° a> % (14) 


If rp is assumed to be large, making * negligible 


in comparison to —1, Eq. 14 becomes 


+ « 6) 


rh... 


Y 
p? -Py* = + Imp + 
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An analysis of Eq. 10 shows that it is possible to 
use the quadratic equation without resorting to 
computer solutions for determination of C values. 
C values determined from solutions to unsteady- 
state Darcy flow equations can be used for quadratic 
flow due to the fact that the quadratic term in the 
quadratic flow equation is of importance only in the 
region near the wellbore. For the moment let us 
consider W as the mass flow rate at any radial 
location in the reservoir. Since W is at its maximum 
value near the wellbore and since cross-sectional 
area perpendicular to flow (A) varies with the 

Bw? 
4n272h2 
negligible except near the wellbore. It follows, then, 
that a pressure profile for quadratic flow in the portion 
of a reservoir away from the immediate vicinity of 
the wellbore is nearly identical to that for Darcy 
flow at the same flow rate, even when there is 
considerable deviation from Darcy flow near the 
wellbore. Since C is determined from material- 
balance considerations, since the region near the 
wellbore contributes a negligible amount of mass to 
the total flow after the first few seconds of reservoir 
life and since deviation from darcy flow occurs only 
near the wellbore, C values determined from studies 
of Darcy flow can be used for quadratic flow calcu- 
lations. 


distance from the wellbore, the term is 


Calculation Procedures 


For the early transient period Eq. 15 becomes 





lnr ~ ra a aie ee 


As before, paye occurs at Tgyep, which is again 
defined as Crp so that 


Pave? ie Pat -7 
= InCryp = 4 mi) 





lnrgueD 


Eq. 5 again applies for material-balance purposes. 
For the post early transient period, 


ae 
Pave Py if . (18) 





ln guedD = InCrep = Y/L 


The calculation procedures for quadratic flow are 
entirely analogous to those for Darcy flow. One 
item of theoretical importance arises from neglecting 
the term 1/rp in Eq. 14. For small values of rp the 
term 1/rp can have an appreciable value. In such 
a case, the use of Eq. 15 would not be justified. For 
most practical situations, however, the use of Eq. 
15 will not cause appreciable error because rgp 
will be large for any time period of practical interest. 
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SUPERPOSITION 


Data on flowing gas wells usually are reported 
in the form of back-pressure tests. For tests of 
short duration, the back-pressure test consists 
of single-point readings of well pressures for super- 
imposed rates. Therefore, any effort to reproduce 
test data must include the effects of superposition 
in the calculation procedure. 

Only the shut-in well pressure, top- and bottom- 
hole pressures at the end of a flow period, flow rates, 
gas gravity, reservoir temperature, formation thick- 
ness and duration of flow are reported for back- 
pressure tests. It is necessary, therefore, to deter- 
mine values for k, 8B, @ and r,, by trial and error. 
When the trial values allow reproduction of the test 
results, it will be assumed that they are correct. It 
is believed that the values will be relatively correct 
if, and only if, the number of data points equals or 
exceeds the number of constants to be determined. 

Superposition effects are calculated in the follow- 
ing manner. 

1. In their order of occurrence, label each of the 
rates imposed upon the well by a number. The first 
rate is Rate 1, the second is Rate 2, etc. 

2. Calculate the \ rates which compose a total 
rate n by the equations, 


A rate n=raten-—raten-l ....... (19) 


and 


rate n = MORER «2 tee se» « OD 
Thus, Rate 1 = A Rate 1, Rate 2 = A Rate 1+ A 
Rate 2 = Rate 1 + Rate 2 — Rate 1, etc. 

3. Calculate the production associated with each 
A rate n at any time ¢ by multiplying A rate n by 
the duration of time during which A rate n has been 
in effect. 

4. Using the Darcy flow equations, calculate p,? 
— Py? for each A rate at every data point. A data 
point is defined as a point in time where rate, 
duration of rate and p,, are known. 

5. For each rate, calculate the effect of the 


quadratic term # in Eq. 15. 


6. Add the individual p;2 — p,,2 values to the 
quadratic effect to obtain the total pressure drop 
at any data point. 


RESERVOIR EXTENT TEST 


An analysis of the C values shown in Fig. 1 
suggests that a method for determining the extent 
of a reservoir during the early transient period can 
be developed. ryp and C can be determined by 
using Eqs. 16, 17 and 5. If C is greater than the 
minimum value of .409 (it is no longer constant), 
thenrgp/Tep can be read from Fig. 1. Mathematically, 
tap /t ep ={ where f is the value read from Fig. 1, and 


“4D Lon ae 


If C = .409, f can have any value from 0 to approx- 
imately .5. To determine the extent of a reservoir, 
therefore, it is necessary to run a test at least until 
C is greater than .409. In actual practice it would 
be advisable to run the test until C was considerably 
above .409 to make sure that the change was not 
due to the random deviations expected in any at- 
tempted calculation of reservoir performance, 

For quadratic flow, the steps in performing this 
calculation are as follows. 

1. Calculate rgp using Eq. 16. 

2. Using rgp from Step 1 and G,, from test data, 
calculate Apgye using Eq. 5, and Paye from Apaye 
= Pi — Pave 


3. From Eq. 17, 


esis rs i" -L 





InC + 1 = 
nC + In rgp "773 
Calculate C using this equation, pgye from Step 
2 and rgp from Step 1. 

4. From Fig. 1 read /, typ /rep for C calculated 
in Step 3. Then calculate rep using Eq. 21. 


DISCUSSION 


To use the new calculation method, the engineer 
must know permeability, porosity, well radius, 
reservoir thickness, gas properties, outer drainage 
perimeter and f. Although permeability, porosity, 
well radius and B could be determined by other 
means, it is felt that better results will be obtained 
if they are calculated from well-test data. 

The well tests can be run at any time during a 
reservoir’s life if the reservoir is previously allowed 
to build up to a static pressure. However, one of 
the primary advantages of the new method is that 
it allows predictions of reservoir performance be- 
fore the reservoir pressure has declined; therefore, 
to gain the fullest benefits from well tests, they 
should be made as early as possible. 

After the engineer has obtained the necessary 
data, the new method provides for him a way to 
calculate the total production associated with any 
well pressure when the well is being produced at 
a constant rate. This is accomplished without iter- 
ation and without calculating production for any 
well pressure other than the one of interest. 

It is realized that the necessity of determining 
k, B, & and r,, by trial and error limits the useful- 
ness of the new method at present. However, develop- 
ment of digital computer programs to perform the 
lengthy trial-and-error calculations should result 
in use of the new method for routine field calcu 
lations. 

No rules have been given for determining the 
averaged values of z and pto be used in a given 
calculation. It is felt that zg,- and p should corre- 
spond with Pgye, but the use of values correspond- 
ing with p; during the ‘‘early transient period” or 
Pe during the “post early transient period’’ probably 
will not cause appreciable error. 
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The analysis justifying the use of C values 
determined from Darcy flow studies for quadratic 
flow calculations should also apply for z and p 
when they are variable functions of pressure. To 
use C values determined from Darcy flow studies 
for another type of flow, it is necessary that the 
pressure profile in the reservoir away from the 
immediate vicinity of the wellbore for the flow in 
question be almost identical to that for Darcy flow. 
This is true for variable z and p because most of 
the pressure drop, and thus most of the deviation 
from Darcy flow, occurs near the wellbore. 


NOMENCLATURE 
A = area cross-sectional to flow 
C = rave/Td =TaveD/aD 
f = Taped 
G p = dimensionless cumulative production 
G,, = Cumulative production, mol 
G,, total = total number of moles in the reservoir 
Gp = cumulative production, Mscf 


h = net reservoir thickness 


Bw? 
~ 2r,, b2 
a 
k = permeability 
[ae = 2 
Zaye RT 
M = gram molecular weight 


N = flow rate in gmmol/sec/cm of reservoir 


thickness 
p = pressure 
Pave = avetage pressure 
Pp = dimensionless pressure 
Pe = pressure at the outer drainage radius 
py = formation pressure 
p; = initial pressure 
)jp = dimensionless initial pressure = 1 
Py = well pressure 


Pwp = dimensionless well pressure 
\Pave = average pressure drop 
\Pavep = dimensionless average pressure drop 
q = volumetric flow rate 
R = gas constant 
r = radius 
Tave = tadius at which Paye occurs 


TaveD = dimensionless radius at which Paye 


occurs 
rp = dimensionless radius 
tq = drainage radius 
Tap = dimensionless drainage radius 
To = outer drainage radius 


‘ep = dimensionless outer drainage radius 
well radius 


~ 
" 


dimensionless well radius = 1 


TwD 
T = absolute temperature 
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¢ = time 
tp = dimensionless time 
u=Inrp 

W = mass flow rate 


Wp = dimensionless mass flow rate 


Wp = average dimensionless mass flow rate 
Y = Wyu/kh 
z = compressibility factor 

Zave = average compressibility factor 

y, @ = back-pressure-test constants 

B = unnamed empirical constant 

= porosity 

p = density 


pL = viscosity 
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APPENDIX 


DIMENSIONLESS VARIABLES 


The computer solutions used for determining C 
values are provided in dimensionless form. The 
dimensionless terms are dimensionless pressure, 
dimensionless rate, dimensionless time 
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where G,, total is total moles of gas initially present 
in the reservoir. From Eq. 2, 

















(1 - *) 
rap = exp Pwd sis i ee ci (24) 
ZaveWp 
and 
2 (l« 2) 
a, = exp cep » «i 
ZaveWpD 







Total number of moles is as follows. 











P;o7 (7.2 ~ Ty 2)b 


G,, Total = 
” ZayeRT 





. (26) 












If r,,2 is assumed negligible in this equation (this 
is valid for almost all practical problems involving 
well behavior) and in Eq. 5, 








Gp = G,,/G,, Total 
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where Apguep is dimensionless average pressure 
drop. Also, since 


Apave = PP; — Pave 


then 


AbaveD = 1 — Pavep- ee ae eer er ae 


Thus, 





z 
TdD ; 
Gp =(1- Paved) yp + + + + ee + (30) 
TeD 
and 
Gpren2 
PaveD = 1 ee » + & © . , ma eS OS (31) 
rap” 





From Eq. 6, 


a Pavep*? — 1 . 


ZaveWp 


Substitution of Eq. 31 into this gives 
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bs Sorep) 
: . 
TdD 
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Zave"p 





Further substitution of Eq. 25 into Eq. 32 produces 





2 
1 Gnr.n 2 
a ., D*‘eD 5 ae 
Zave4n 2(1 - Pwd ) 
Se 
Zave"p 


or 





For tie post early transient period, in dimension- 


less terms, 





1 
InCrep = Zave¥p (Paved? = Pwd?) 
and 
2 2 
\ Gprep Gplep 
MaveD = rap? = gh = Gp, 
e 
or 
Paved = 1-Gp. +... toy eR, 28 « wp) 
Therefore, 
1 
InC = [a ~Gpy’- Pwo | ~ lnrep, 
Zavep 
. ££ 2-4 “S@ 22 ee eS we . . . (36) 
or 
1 2 2 
c - exp | : [a - @) — Pwo ]- treo} 
Zave"D } 


eee see owe eae 
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CALCULATIONS 
DARCY FLOW — CONSTANT-RATE CASE 


To illustrate the use of the calculation method for 
a constant-rate case, let k = .3 darcies, rep = 1,000, 
ry = 15.25 cm, p; = 136 atm, Zgye = 1, p = .016 cp, 
T = 332°K, M = 20.86 gm/gm mol, ¢ = .20, R 82.06 
cc atm/gm mol °K, b = 304.8 cm, and W = ~—5848.1 
gm/sec = 20,300 Mscf/D. 

From Eq. 2, 


(n) (.3) (20.86) (136)? (304.8) (,_, 2 
(5858.1) (.016) (82.06) (332) ‘1 ~Pwd ) 





In dD = 


43.478 (1 -p, 7); 


it 


from Eq. 4, 


In raved = 43.478 (Paved? — Pwo); 
and from Eq. 5, 
(1) (82.06) (332) G,, 





\Dave = 
ave (304.8) (.2) (m)(rg2 —r 4,2) 
G 
= 142.257—= . 
rae 
Also, 


Gm = Apavetg 2/ 142.257. 


For the first point let p,, = 126.5 atm. Table 2 
shows calculations using this value. The dimension- 
less drainage radius is 351.4. Therefore, rgp/rep = 
.3514. C is then read from the curve in Fig. 1 as 
409. ravep =-409 rgp = 143.72, and Intgyep = 4.968. 
Dimensionless average pressure is calculated as 
.98967, dimensionless average pressure drop as 
-01033 and average pressure drop as 1.40488 atm. 
rq = 15.25,rgp = 5358.85, and ryp? = 28.7173 x 10°. 
G,, is then calculated as 283,601 gm mol production. 
Time ¢t is then calculated by the equation 


t= oe ee ae a 





1011.6 seconds, or 16.86 minutes. 


For the next point let p, = 125.0 atm. Calculated 
by the same process used for Point 1, G,, = 1,536,110 
mol and t = 5479.2 seconds. 

For the next point let p,, = 120.0 atm. Since largp 
is gteater than In 1,000, the reservoir is in the post 
early transient stage of production. For this situation 
Eqs. 7 and‘8 are used for calculation purposes. G,, 
is then calculated as 9,071,192, and t as 32,356.1. 


DARCY FLOW — CONSTANT-PRESSURE CASE 


The calculation method can be used to calculate 
a constant-well-pressure case (pressure at the sand 
face remains constant) by treating it as a series of 
constant-rate cases. It is assumed that the pressure 
profile at a given rate is accurately described by 
C values determined from constant-rate studies. The 
procedure consists of assuming a series of decreas- 
ing fates, calculating the production associated 
with each rate and dividing the difference in produc- 
tion between two successive rate values by an average 
of these two rates to obtain the elapsed time between 
the occurrence of the two rates. Since this is essen- 
tially a numerical integration procedure, errors will 
occur. The magnitude of the error should depend upon 
the size of the time step and should disappear as the 
number of time steps used approaches infinity. 

The constant-pressure case presented here is 
calculated in dimensionless terms. This is done 
so that it can be conveniently compared with digital 
computer results for the same problem. 

The basic equation for use in the early transient 
period is Eq. 33, which can be written as 


2(1—p ,p) 
1+2ZgyeWp InC =  ~ Corea Sa . 


Zave 
From this is obtained 


2(1 — pp?) 
2 a 
- cod fo Fo VE Faye 


Zave 











and 
G 1 — 1+ ZapveWpInC (39) 
D= oe «6 aecdie 
2 2(1-p,,p?) 
: fi em 
Zave Wp 


The corresponding equation for the post early tran- 
































TABLE 2 

Point Pw PwD Pwo 1 ‘ Pwo Inran ‘4D raD/"eD c ‘aveD 
1 126.5 930147 865173 - 134827 5.8620 351.4 23514 -409 143.72 
2 125.0 291912 844782 -155218 6.7486 852.8 -8528 .4395 374.81 
3 120.0 -88235 -778542 -221458 9.6286* 1000 1.000 .470 470.00 

\ 

Point Inr, veD PaveD PaveD APavep APave ‘a rg/ 10° G,, t 
1 4.9680 97944 -98967 -01033 1.40488 5358.85 28.7173 283601 1011.6 
2 5.9264 -98109 -99050 .00950 1.2920 13005.2 169.1352 1536110 5479.2 
3 6.1528 292006 -95920 .04080 5.5488 15250.0 232.5625 9071192 32356. 1 


*Since this value is larger than In 1,000, the reservoir is in the post early transient stage of production. 
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sient period is Eq. 36, which can be arranged as 


(InC + lnrep) ZeueWD + Pwo =(1 — Gp)? 
which yields 





Gp =l- ¥zaveWD (InC + lnrep) + non” e . (40) 


For calculation purposes, let p,,p =.1, rep = 1,000, 
lnrep = 6.908, and zgye =1.0. 

During the early transient portion of the calcula- 
tion, tgp /T ep can be calculated by the equation 


1 2 
J ne —PwD 
TdD’"eD = ©XP pe Wp TeD 
ave 


» « « (41) 


(see Eq. 24), and C can be read from Fig. 1. The 
dimensionless time interval Atp associated with a 
production interval is calculated by the equation 


AGp(r ep 2 = T wD?) _ AGp x 106 
W 1 





D "p 


where 


Wp at start of interval+Wp at end of interval 





(43) 


To calculate the first time step, it is necessary 
to assume an average rate. If the first rate chosen 
is high enough, this assumption produces a small 
and constant error which is of importance for only 
the first few time steps. 

Tables 3 and 4 show early transient and post 
early transient calculations, respectively, for the 
constant-well-pressure case. Fig. 2 is a graph of 
these values of rate and production vs time, along 
with values from a digital computer solution to the 
same problem. Except for tp values of less than 100, 
agreement seems quite good, especially when the 
approximate nature of both calculations is considered. 

Values of tp less than 100 correspond to small 


tap values and, consequently, to very small values 
of time. As explained earlier, there is a loss of 
accuracy at small values of rgp. For this reason 





TABLE 3 — CONSTANT-WELL-PRESSURE CASE, EARLY TRANSIENT 


(2) (3) 
(1) (-Pyp*) rap (4) = (5) 


a Wp InC 


(6) (7) (8) 
1+WpInC 1-(6) Atnrep-(2)] 


Wp Wp 'eD c 





(9) (10) (11) (12) (13) (14) 
e2'8) Gp AGp Wp Atp/10°® th 





80 1.238 
70 1.414 
-60 1.650 
50 1.980 
-40 2.475 
035 2.829 
30 3.300 
28 3.536 
26 3.808 
224 4.125 
022 4.50 

20 4.95 

18 5.50 


-00342 .409 -.715 .534 
000412 .409 -.626 .611 
-00522 .409 -.536 .680 
-00725 .409 -.447 .743 
0119 + .409 -.358 .801 
0169 .409 -.313 .828 
-0270 .409 -.268 .855 
034 .409 -.250 .866 
045 .409 -.232 .876 
062 .409 -.215 .886 
090 .409 -.197 .896 
-142 =+.409 -.179 =.905 
245 .409 -.161 915 
16 6.188 .485 .409 -.143 .926 
015 6.60 0735 4 =.422 -.1296 .9329 
-1433 6.908 1.0 -470 -.1082 .9443 


464 
+389 
+320 
+257 
199 
0172 
145 
134 
124 
114 
104 
095 
-085 
-074 
-0671 
0557 0 


11.340 
10.988 
10.516 
9.856 
8.866 
8.158 
7.216 
6.744 
6.200 
5.566 
4.816 
3.916 
2.816 
1.440 
616 


*Assumed. 


-553x 10-® .553x 1075 1.0* 
5.90104 .659x107> 1106x1075 .75 
3.68104 .869x107® .210x10°°  .65 
1.902 104 .135x 1074 .482x1075 .55 
7.07x 102 .281x107* 1146x1074 .45 
3.48x10? .494x1074 .213x107* .375 
1.355 10° 107x107? .576x1074 .325 
8.47x102 .158x107> .51x 1074 ~=.29 
492 6252x107? .94 x 1074 
261 437x107? .185x 1073 
123.1  .845x107? .408x 1075 
50.2 .189x107? .1045x107? 
16.7. 509x107? .320x 1072 
4.22 .177x107' .126x 107! 
1.85  .363x107' .186107' 
1 557x107" .194x107' 


8.39 x 104 .553x 1075 5.53 
.141« 1075 6.94 
.323x 107° 10.2 
876 x 1075 18.9 
324 x 1074 51.3 
568 «1074 108 
.1772 1073 285 
1759 x 1073 461 
.3481x 1073 809 
.740 «1073 1550 
.1779x1072 3320 
.4976x1072 8300 
.1684x107' 25140 
.7418x107' 99320 
.1198 219100 
1325 351600 





POST EARLY TRANSIENT 


(2)* 
(1) ZayveWp(In C + Inrp) 
Wp + PwD 
14 -87142 
13 -80989 
012 -74836 
1} -68683 
-10 -62530 
.09 «56377 
.08 -50724 
.07 .44071 
06 37918 
05 031765 
.04 225612 
-03 219459 
02 213306 
01 -07153 





-7908 
-7508 
7087 
+6639 
-6158 
+5636 
-5061 
4411 
+3648 
+2675 


(5) (8) 
AGp 'p 


-0108 427800 
-0335 675900 
0349 955100 
-0363 127 1000 
-0380 1633000 
-0400 2054000 
0421 2549000 
+0448 3146000 
-0481 3886000 
+0522 4835000 
+0575 6113000 
-0650 7970000 
-0763 11020000 
0973 17510000 


*C = .470; In C = -.755; In rep = 6.908; and In C + Inrgp = 6.153. 
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no values were computed for Wp values greater 
than .80. Fortunately, the production at this point 
is so small that higher rates can be neglected with 
little effect on the over-all calculation. This neces- 
sarily will be true in all cases involving a constant 
well pressure. 


QUADRATIC FLOW — CONSTANT-RATE CASE 


To illustrate a quadratic flow calculation, the same 
case that was presented for Darcy flow, with the 
addition of B = ~ .648 atm sec2/gm, is recalculated 
using quadratic flow equations. 


From Eq. 12, 
vit « WZ aveRT 
khnM 





__+ (5848.1) (.061) (82.06) (1) (332) 
(.3) (304.8) (7) (20.86) 





























= 425.414, 
and 
BW2z0,eRT 
J/L nibs... tae 
2a 2, b2M 
+(.648) (5848.1)? (1) (82.06) (332) 
2(m)? (15.25) (304.8)? (20.86) 
= 1035.0. 
és CONSTANT PRESSURE CASE Pyo=.! reo 1000 
CALCULATED Wy a 
© DIGITAL Wp 
0.5 tT] 6 CALCULATED Gp | om 
> DIGITAL Gp " 
0°04 > + } x 
°o b J 
Fy 0.3 © 5 
. © 7 : 
z - oe - 
0.2 Or $ 4 
© Bt " 
" Bop 
0.1 + + + + a Fo + 
| 6 é “ta 
0.0o-5—s—A—s-#-a-n—an—oi__6 4 Q 
10 102 103 10% 10° 108 107 10° 


FIG, 2 — COMPARISON OF HAND-CALCULATED 

VALUES OF 4 AND Gp WITH COMPUTER VALUES — 

DARCY FLOW (CONSTANT-PRESSURE CASE pywp = 
ys %eD = 1,000), 








— ANALOG 








4 CALCULATED 





PRESSURE (ATM) 


8 
| 
| 




















4 ba 
° 80,000 160,000 240000 320000 400000 480,000 
TIME (SECONDS) 


FIG. 3 — COMPARISON OF ANALOG AND CALCULATED 
VALUES OF fw» VS TIME — QUADRATIC FLOW. 


DECEMBER, 1961 


Eq. 16 then becomes 


p;? — py? — 1035.0 
aD ~ 425.414 





As before, Eq. 5 becomes 
Abave = 142.257 G,,/rq? 
or 


Gm = Npavyetd’/ 142.257; 


and, from Eq. 38, 
t = G,, (20.86) / 5848.1. 


Fig. 3 shows a comparison of calculated values 
of py vs t for this case with values computed on an 
analog.10 The agreement is very good considering 
the problems associated with analog computations. 
The observable deviation in values at very small 
values of time is due to the inability of the x-y 
plotter to follow the analog computations at the 
particular time scale used to obtain this figure. 
The deviation at large values of time is probably 
due to analog drift. 

Fig. 4 shows a comparison of the calculated 
values of p, vs ¢t for Darcy and quadratic flow 
given in Tables 2 and 5, respectively. As would be 
expected, a lower well pressure is indicated at a 
given time for quadratic flow. 


QUADRATIC FLOW — CONSTANT-PRESSURE CASE 


A constant-pressure case involving quadratic 
flow can be calculated in a manner analogous to 
that used for Darcy flow. The basic equations to be 
used here are Eqs. 16,17 (or 18), and 5. After AG,,, 
the production during an interval, is determined, the 
elapsed time during the interval is calculated by the 
equation 


MAG,, 





At = (44) 


W 


132 4 — DARCY 
© — QUADRATIC 






s 
@. 


PRESSURE (ATM) 
x) rx) 
° > 


12 
4,000 12,000 20,000 28,000 36,000 


TIME (SECONDS) 
FIG. 4 — COMPARISON OF CALCULATED VALUES OF 
bw VS TIME FOR DARCY AND QUADRATIC FLOW. 
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where 


y= W_at start of interval + W at end of interval 
2 





~ » (45) 


. 
7 
. 
. 
. 


SUPERPOSITION 


An actual field case is used to illustrate applica- 
tion of the method when superimposed flow rates 
are involved. The data used here are those of Baumel 
and Breitung.® These data also were used by Cornell 
and Katz? in a paper in which they developed a 
gtaphical calculation procedure and attempted to 
correlate the data. 

As a first trial, let 6 = .20375, k = .10125 darcies, 
B = — 2.7652 atm sec2/gm, and r,, = .22969 ft = 
7.0009 cm. 

These values were obtained by a computer method 
which is not described here. It is recommended that 
the engineer does not try to perform a trial-and-error 
procedure to determine ¢, k, 8 and r,, unless he has 
access to a digital computer. 

The test data include: 4 = 21 ft; p; = 3,127 psia 
(24-hour shut-in pressure) = 212.721 atm; gq, = 3,710 
Mscf/D for 180 minutes, p,, = 3,087 psia; q2 = 5,980 
Mscf/D for 200 minutes, b,, = 3,059 psia; 93 = 8,191 
Mscf/D for 180 minutes, p, = 3,035 psia; q4 = 
14,290 Mscf/D for 190 minutes, p,, = 2,942 psia; T= 
370°K; M = 20.3 (gas gravity = 0.700); p = 0.021 
cp (average value); zg,,- = 0.87; and casing diameter 
= 5.5 in. 

The equations to be used are Eas. 2, 3, 4, 5, 6, 
16, 19 and 20. In Eq. 16, Y/L and J/L will be 
Have RT 2 BzayeRTM 
a Bel afl” Kem 

ak 272 r,, 


respectively. N is negative and equals W/MhA, which 
is rate in gram moles per secend per centimeter of 
reservoir thickness. Therefore, 


rewritten as —N 








Tep is unknown, but for a 12-hour test it is doubt- 
ful if rgp/tep would much exceed .50. If rgp/rep 
is less than .5, C is constant. Therefore, it will 
be assumed that C and InC have constant values of 
410 and —.8915, respectively. 

From Eq. 6, 


Pave? = p;2 + (Y/L) InC. 


Since it has been assumed that C is a constant, 
Pave is constant at a constant flow rate. It follows, 
therefore, that there is a constant Ap,,,, associated 
with each A rate. 

Before proceeding, it is necessary to calculate 
rates, A rates, production associated with the A 
rates, values of J/L and Y/L associated with the 
rates and A rates, respectively, and Apagye associated 
with each A rate. Now recall that 


Seuq GghkT 
hd 7 rap Ty? 





\Pave sa 


and that 


In tfavep = InCrgp 


—7k 2 2 
= noes ee oe 
Nu ZeygRT (Pave Pu 





3 Riee* oe ty * 
Y/L 


Also, 





| zendeT 
dD ~ 
bd 7 \Pave ’ 


















































so that 
—N (.021) (.87) (82.06) (370) . 
a a ae 7 : N = 
“Vho 7 \pave lu" — ail 
and or 
—2.7652) (.87) (82.06 70) (20. ‘ west 
T/L =n? {2:7652) (.87) (82.06) (370) (20.3) 6 cuss suing a 
2(3.1416)? (7.0009) Py = Pave ~ ~ in 2 
2 hd m Ntave Tw 
= 10,729.8N ¢ 
TABLE 5 — QUADRATIC FLOW CALCULATIONS 
Point Ps Pe P;7~ — p?- 1035-p,? Inrup rap 'ap/"eD c Cased 
1 122.5 15006.3 3489.7 2454.7 5.7702 320.6 +3206 -409 131.13 
2 121.0 14641.0 3855.0 2820.0 6.6288 756.5 7565 +424 320.76 
3 116.0 13456.0 5040.0 4005.0 9.41*** 1000 1.0000 -470 470.00 
Point InrayveD Pave Pave APave ‘d rf /\0° GC, t 
1 4.87615 18115.7 134.595** 1.405 4889.15 23.90379 236086 842.1 
2 5.77073 18131.0 134.652 1,348 11536.63 133.0938 1261168 4498.5 
3 6.1528 17108.5 130.800 5.200 15250.0 232.5625 8500976 30322.1 


*p? = 136° = 18,496. 
**This is identical to the value for Darcy flow with p,, = 126.5. 
***Greater than In 1,000; therefore, post early transient stage. 
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2D 
ill 


reD 


1.13 
).76 
).00 


2.1 
18.5 
12.1 








Data Rate 

Point (Mscf/D) N 
1 3710 ~ .08013 
2 5980 -.12916 
3 8191 -.17691 
4 14290 - .30864 


Data Points (G,) 





TABLE 6 — RATES, A RATES AND ASSOCIATED VALUES 








AN 1 2 3 4 
- 08013 464 979 1443 1932 
- .04903 0 315 599 899 
-.04775 0 0 277 569 
-.13173 0 0 0 804 


*p,*= 45,250.2. 











AN N? J/L 

- .08013 -0064208 68.894 
- .04903 -0 166823 178.998 
~- 04775 0312971 335.812 
-.13173 -0952586 1022. 106 

(Y/L) InC, or — AP ave 

¥7& Pave-~ p,*" Paus (atm) (atm) 
139.739 - 124.6 45125.6 212.428 2293 
85.503 - 76.2 45174.0 212.542 179 
83.271 - 74.2 45176.0 212.546 175 
229.724 - 204.8 45045.4 212.239 -482 





Since 


Gp (453.6) (103) 
si 380.15 





where Gp is production in thousands of standard 
cubic feet, the equation becomes 


ZaveGPRTC7(453.6 x 10%) 
bd 7 AP ave 2 (380.15) 





PZ = Days—1/2 Y/L In 


= Pave? - 1/2 Y/L 


.87 Gp (82.06) (370) (.41) 7 (453.6 x 10°) 





6A0 (.20375)(3.14.159)(Apave) (7.0009)? (380. 15) 
Gp 
MPave’ 





= Pave - 1/2 Y/L In 263.878 


Thus, from Darcy flow equations, pi - aK, terms 
associated with each A rate at each data point are 
calculated by the equation 


2 
Ap ave 


The sum of these and the J/L term at a point gives 
the total p;“ —p, for the point, from which can be 
obtained the well pressure. 

Calculations leading to the determination of p,, 
values are presented in Tables 6 and 7. Considering 
all factors, agreement between calculated values of 
Py and field values is quite good. Maximum deviation 
of a calculated value from the reported value is 
3.4 psi for the third data point, which is well within 
experimental error. The third data point proved 
difficult to fit. If the pressure at this point was 
lower than reported, or the rate higher, an excellent 
fit could be obtained. 

Difficulties in matching reported data might occur 
in cases where reported pressures and rates are 
inaccurate, where the formation has been damaged 
with accompanying loss of permeability, when the 
time consumed in changing chokes is great enough 
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oF = py? op? ~Peve? + 1/2 Y/L 1n 263.878 ~e— 


to allow considerable pressure change, or if major 
variations from radial geometry are present. 


RESERVOIR EXTENT TEST 


For this example assume that the reservoir treated 
in the quadratic flow constant-rate case is being 
tested. In the present case the production data are 
known and rep is the unknown which is to be deter- 
mined. Let p,, = 121 atm, G,, = 1,261,168 gm mol, 
W = -—5,848.1 gm/sec, and t = 4,489.5 seconds. 

This is a reverse calculation of the second point 
in Table 5. Y/L again equals 425.414 and J/L = 
1,035.0. From Eq. 16, 


2 — Py? — 1035.0 





























lnrgp = 
ad 425.414 
855.0 — 1035.0 
oa 35:0 _ 6.6288. 
425.414 
TABLE 7 — DETERMINATION OF p,, VALUES 
263.878 
\Rate P?- Paved V/2AY/L) es 
1 124.6 69.870 900.68 
2 76.2 42.752 147.418 
3 74.2 41.636 150.787 
4 204.8 114.862 547.46 
(1) (2) (3) (4) 
Gp Gp(263.878)/\p, ve In(2) V2AY/L)In(2) sp? -p,2 
\Rate 1 
464 417916 12.9431 904.33 1028.9 
979 881766 13.6898 956.51 1081.1 
1443 1299681 14.0776 983.60 1108.2 
1932 1740114 14.3696 1004.00 1128.6 
A Rate 2 
315 464367 13.0485 557.85 634.0 
599 883034 13.6916 585.34 661.5 
899 1325288 14.0972 602.68 678.9 
\Rate 3 
277 417680 12.9425 538.87 613.1 
569 857978 13.6624 568.85 643.0 
ARate 4 
804 440158 12.9950 1492.63 1697.4 
\Rotes and i 
i Bes Vv. 2. 8 
Diese Repganved As¥ociated p*- p,~ Values prop. Py Pw 





Point p,, (psia) 1 2 3 4 J/L Total (atm) (psi) 


1 3087 1028.9 0 0 0 

2 3059 1081.1 634.0 0 0 179.0 1894.1 108.221 3060.8 
3 3035 1108.2 661.5 613.1 0 335.8 2718.6 206.232 3031.6 
4 2942 1128.6 678.9 643.0 1697.4 1022.1 5170.0 200.200 2942.9 


68.9 1097.8 210.125 3088.8 
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rap = €9:6288 = 756.5, 


rq = 15.25 (756.5) = 11,536.63 
and 


rq2/106 = 133.0938. 


From Eq. 5, 


APave = 142.257 5 


Td 
is 142.257 (1261168) 
133093800 





= 1.348 atm, 
and 
Pave = 136 — 1.348 = 134.652 atm. 


From Eq. 17, 


- P.? - J/L : 
Y/L ee 
__18131.0 — 14641.0 ~ 1035.0 
425.414 


InC = Pave : 








= — .8580, 


and 
C = e~-8580 = .4240. 
From Fig. 1, 
tap ‘Tep = -7565, 
and 
Top = %gp/+7905 = 756.5/.7565 = 1,000. 
Also, 


r 


— 6.6288 


© = leDlw = 15.25 (1000) = 15,250 cm = 500 ft. 
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An Approximate Method for Determining Areal Sweep 


Efficiency and Flow Capacity in Formations with 
Anisotropic Permeability 


M. MORTADA 

MEMBER AIME 

G. W. NABOR 

JUNIOR MEMBER AIME 


ABSTRACT 


The effects of anisotropic or directional permea- 
bility on the areal sweep efficiency and the flou 
capacity are examined. The paper points out the 
importance of taking directional permeability into 
consideration in planning a flood. It analyzes the 
two-dimensional flow pattern associated with the 
skewed line drive for a unit mobility ratio. The direct 
and staggered line drives are treated as special 
cases of the skewed line drive. 

Analytical expressions are developed for the 
areal sweep efficiency at breakthrough and the 
flow capacity. They are related to the spacing 
between like wells, the distance between a row of 
injectors and the nearest row of producers, and 
the degree of skewness of the line drive. The 
latter quantity is defined such that it is equal to 
zero for the direct line drive and equals one-half 
for the staggered line drive. The areal sweep 
efficiency and the flow capacity depend also on 
the orientation of the flood pattern with respect to 
the principal axis of anisotropy. 

The paper provides a simple method for determin- 
ing the areal sweep efficiency and the flow capacity 
for a formation in which the permeability in the 
bedding plane is anisotropic. 


INTRODUCTION 


Directional or anisotropic permeability is mani- 
fested by the ability of the formation to conduct 
fluids more readily along certain preferred directions. 
This situation occurs in many producing forma- 
tions!-4 and is usually attributed to depositional 
features in which the sand grains are oriented in a 
preferred direction. In some cases it results from 
the formation of a major and a minor fracture system. 4 

Directional permeability should be taken into 
account in many phases of the production and 
exploration activities. Recognizing its existence 
in the formation of interest and planning accordingly 
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can lead to increased recovery and substantial 
savings. For instance, the areal sweep efficiency in 
a water flood depends to a great extent on the orien- 
tation of the flood pattern with respect to the prin- 
cipal axis of permeability. Anisotropic5 permeability 
is specified by the directions of its three principal 
axes and the permeability along each axis. The 
principal axes of permeability are mutually perpen- 
dicular. 

This paper deals with the areal sweep efficiency 
at breakthrough and the flow capacity for formations 
with anisotropic permeability. The flood pattern 
considered consists of alternate rows of injecting 
and producing wells. The rows of wells are parallel 
and form a developed, skewed line drive which is 
illustrated in Fig. 1. The staggered and direct line 
drives are treated as special cases of the skewed 
line drive. 


PRELIMINARY CONSIDERATIONS 
The formations considered in this paper are 
assumed to be of uniform thickness and porosity. 
The permeability is uniformly anisotropic, with 
each principal axis of permeability maintaining 


y=d 


aan pee, 
| 
ik * whos 
| e 1 ln 
ae es nen 
| ! 
ae 
y=0 e L--e--J- e x axis 
a 
°o oO °o ° 


e INJECTION WELL 
o PRODUCTION WELL 


FIG. 1 —- THE DEVELOPED SKEWED LINE DRIVE. 
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the same direction and magnitude at every point 
in the formation. Two of the principal axes lie in 
the bedding plane, and the third axis is perpen- 
dicular to it. Thus, neglecting gravity forces reduces 
the study to a two-dimensional flow problem. It is 
further assumed that the mobility of the fluids is 
the same ahead and behind the front. 


TWO-DIMENSIONAL FLOW EQUATIONS 


The flow velocities in anisotropic permeability 
systems can be described by Eqs. 1. In these 
equations the major and minor axes of permeability 
(n, 2) serve as the co-ordinate system. 








__ ky 5p ir SV 
‘” ph 8m wa 
k dp SV 
Paes . 
= SBE a i 


ky and k¢ are the permeabilities along the major 
and minor axes of permeability, vp, and v¢ are the 
microscopic flow veldécities along these axes, p 
is the pressure distribution, V is the stream function, 
@ is the porosity and yp is the fluid viscosity. The 
flow velocities in another co-ordinate system (l,m) 
which forms an angle 6 with (7,¢) are given by Eqs. 
2. 


Ky cos’ 8 +k; sin’@ Sp 














ail ud 51 
kr sin@cos@-k, sin@cos@ 8p 
Boe 8m 

se 4 2 

, k, sin Atk cos@ 8p 

saints u > 8m 
Ky sin 8cos 8 -k, sin 8 cos 6 Sp 

_ re 4, 

Be 8! 


Eqs. 2 show that isotropic permeability is a special 
case in which the permeabilities along the principal 
axes are equal. For isotropic permeability, Eqs. 2 


k 


k 
become v7 = “ud I? Um= ud P regardless of 


the choice of the co-ordinate system. 

The flow of incompressible liquids in a formation 
that is uniformly anisotropic is described by the 
differential equations: 

2 2 
eS ee ee 





wp 87° =p BE? 
k, SV kp BY. 
Y* ‘Sis 
pp 8m pp 8C 
SF k, 8p 8F 
ss ke 
Ke 8p 8F 


“ae at ar ** - bee 





F (n,¢,t) = 0 gives the position of the flood front as 
a function of time.® 

For the explicit solution of these differential 
equations, it is necessary to specify the pressure 
or the influx (stream function) on the boundary of 
the region of interest, the initial position of the 
flood front and the location of injection and produc- 
tion wells. The solution of these equations is 
simplified by introducing the transformation® xX = 
1/Vhy y= ke and applying it to the differential 


equations and the boundary conditions. The resulting 





equations, 
~ oe 
‘ Bp 8K dy 
~ +e oe, 
. wp BY 8x 
2 2 
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8x sy 
De ae ae 
-.' oar * 
8x sy 
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St pp 8K 8X 
| 8 SF 
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describe the pressure distribution, the stream func- 
tion and the flood advance in an isotropic system 
with permeability equal to Vy ke. The porosity 
and fluid viscosity are equal to those in the original 
system, and the geometry is related to the geometry 
of the original system by the transformation cited 
previously. The solution of Eqs. 4, accomplished 
either by analytical methods or by using a potenti- 
ometric analyzer, yields the areal sweep efficiency 
and the flow capacity for the new isotropic and 
the original anisotropic systems. This results be- 
cause the areal sweep efficiency and the flow 
Capacity are preserved in the transformation. 


ANALYSIS OF THE FLOOD 
PATTERN — AND RESULTS 


In analyzing the flood pattern it is useful to 
examine the pressure distribution due to an infinite 
atray of point sources of strength qg, with spacing a 
and forming an angle 6 with the major axis of per- 
meability. The infinite array upon transformation 
remains a straight line, but with spacing 


—— (28 sin? @ 
aq-=a + 
K» ke 


and forms an angle 


Q = ton’ (./k, 7k tan @) 








SOCIETY OF PETROLEUM ENGINEERS JOURNAL 














with the x axis. The pressure distribution in the 
new system is given by® 


— 
‘anh, /kn Ke 


cosh (22) (9 cos@-X sin 6) 


‘In : a 
—cos cale cos §+f sind) 


Eq. 5 points out two important features. First, at 
a distance from the array >0.6da, the isobars can be 
approximated by straight lines parallel to the array. 
This is because the cosine term on the right-hand 
side of Eq. 5 becomes negligible with respect to 
the hyberbolic cosine. Also, at a distance > 0.64, 
the pressure gradient becomes normal to the array 
and approximately equals 


(ua/28n/k ke) 


Second, the isobars are essentially concentric 
circles inside the circular regions with center at 
the point sources and radius < 0.05a. 


AREAL SWEEP EFFICIENCY 


To complete the flood pattern, let dbe the distance 
in the original anisotropic system between one 
array of injectors and the closest array of producers. 
The corresponding distance in the transformed 
system is 


d= WAG cos @ +k, sin’ @ 


and 


(4/3) = (4/o) | Haake H the 6056+ kysin?a) 





For the same spacing (d/a) inthe anisotropic system, 
(d/a) increases as @ decreases and reaches a max- 
imum of 


(/ a) = (4/0) Wk, Ze 


at @=0. Therefore, it usually is advantageous to 
arrange the arrays of like wells parallel to the major 
permeability axis. This arrangement results in 
maximum (d/a) and, hence, in larger areal sweep 
efficiency. In the remainder of this paper, it will 
be assumed that one array of wells in the transformed 
system either falls on the ¥ axis or that the co- 
ordinates in the transformed system have been 
rotated on angle 6 to accomplish this purpose. This 
assumption causes no loss of generality; it only 


DECEMBER, 1961 


simplifies the analysis. Let x and y refer to the 
co-ordinate system after rotation. The pressure 
distribution, in the infinite strip bounded by the 
lines y = 0 and y = d/2 in Fig. 1, is equal to that 
due to the x-axis array plus the contributions of 
each pair of arrays equidistant from the x-axis. For 
(d/@) > 1.2, the contribution of each of these pairs 
to the pressure distribution in the strip of interest 
is for all practical purposes independent of x. Thus,® 


“ qu 
Puy) 
ath Ky ke 


2 © 
In [cosh 2 — cos ore pen” 





afcosn 2x 5] [eggy 2xtm dea 
Q a 








in - (6) 
(Armd ) 
— 3 
Further, since y < 0.5 d, Eq. 6 becomes 
Oqry 
Pp «& — ey cosh he cos] . (7) 
hy aah /ky ke . . 


Similarly, the potential distribution in the infinite 
strip bounded by the lines y = d/2 and y =d is given 
by 





q 
Py ” — In [cosh 
Oh Ky ke 


(e/a) is a measure of the skewness of the trans- 
formed flood pattern, and equals zero for the direct 
line drive and 0.5 for the staggered line drive. Note 
that (d/a) > 1.2 is not a severe restriction and is 
usually fulfilled for any pronounced degree of ani- 
sotropy [( /ke) > 5] when the rows of wells lie in 
the direction of the major axis of permeability. 

Using the pressure equation given, the stream 
function and time of travel along a streamline may 
be found. From these, the areal sweep efficiency 
at breakthrough is found to be 





c, cost (r/2)(a/a) Se 


= | : 
(7/2) (d/0) in| cos (3/2) (€/G) 


Details of the analysis are derived in Appendices 
A and B. Eq. 9 provides a simple expression for 
evaluating the areal sweep efficiency provided 
(d/a) > 1.2. In the foregoing analysis the initial 
position of. the front was assumed to be at the 
center of the well instead of an ellipse with major 
and minor axes of Tw/Vk & and Tw / Vey» respec- 
tively. This approximation can introduce some error 
in the areal sweep efficiency and is currently under 
study. 
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FLOW CAPACITY 


In analyzing the flow capacity, the geometry and 
size of the sources and sinks must be taken into 
account. This is true because the pressure gradient 
in the neighborhood of the sources and sinks is 
very large; consequently, an appreciable part of 
the resistance to flow occurs in these regions. 
Therefore, a more detailed examination of the 
pressure distribution in the neighborhood of the 
transformed well is necessary. 

The analysis is greatly simplified by noticing 
that: (1) the isobars are essentially concentric 
circles inside a circular region the center of which 
is the point source and the radius of which is .05 
@; and (2) cylindrical sources or sinks are trans- 
formed into elliptical cylinders, and the isobars 
due to an elliptical source or sink at a distance 
larger than seven times the focal distance / of the 
ellipse are circles for all practical purposes, pro- 
vided the sources do not interfere with each other. 

Hence, the isobars are circular in a portion of the 
transformed reservoir, provided (7/) < (0.05 @). For 
normal well spacings and well radii, this criterion 
is usually satisfied. Choosing a radius 7,, such 
that (7f) <1r,, < (0.05 @) and using this radius as a 
fictitious boundary, the pressure drop is found 
through the main body of the reservoir to 7,,, and 
from 7, to the transformed source or sink. These 


pressure drops are added and the final flow capacity 
equation follows. 


_ 1p Ska ke (Piw™Ppw) (10) 
d i a | _(Inéw 
(oi) elie OF) 
where q is the rate of flow at reservoir conditions, 
and 
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PRACTICAL EXAMPLES 


Shown in Fig. 2 are three proposed patterns for 
water flooding a reservoir which exhibits a large 
degree of permeability anisotropy, 


(ky/ke )=144. 


In these three cases, the rows of wells were oriented 
along the major axis of permeability; i.e., the angle 
6 = zero. Hence, 


(3/0) = (d Johor 


or 


/o) = 12(d/a) , 
and 
(6/4) = (e/a) 


The spacing ratios (d/a) for Patterns A, B and C 
are 4.8, 4.8 and 2.4, respectively, satisfying the 
(d/a) > 1.2 requirement. The patterns were evaluated 
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using Eqs. 9 and 10 and, also, by studying the 
transformed system potentiometrically. The stream- 
lines and flood fronts obtained from the potentiometric 
model were back-transformed to the original reservoir 
co-ordinate system. The maps obtained for Patterns 
A, B and C are given in Figs. 3, 4 and 5, respectively, 
The four fronts shown correspond to times of 0.25 
ty, 0.50 t,, 0.75 t, and 1.00 t,, where t, = initial 
breakthrough time. These fronts, plus the stream- 
lines shown as dashed lines on the pattern, indicate 
the advantage gained by orienting the rows of wells 
along the major permeability axis. The injected 
fluid spreads out rapidly along the major permea- 
bility axis and, thereafter, advances quite uniformly 
toward the row of producers, with high areal sweep 
efficiency resulting. 

The results obtained using Eqs. 9 and 10 are 
presented in Table 1. These results agreed with 


TABLE 1 — RESULTS FOR FLOOD PATTERNS OF FIG, 1 























Pattern 
A a= a oa 
a, miles Vv 1.25 y 1.25 Vv 1.25 
d, miles 0.4y/ 1.25 0.4,/ 1.25 0.2\/ 1.25 
e, miles 0.3\/ 1.25 0.2\/ 1.25 0.4) 1.25 
(d/a) 4.8 4.8 2.4 
(e/a) 0.3 0.2 0.4 
E,, Calculated 0.922 0.915 0.872 
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FIG. 2— PROPOSED FLOODING PATTERNS. 
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the potentiometric results within 2 per cent for 
the flow capacity and within 1 per cent for the 
areal sweep efficiency. The flow-capacity equation 
presented in the table is a reduction of Eq. 10 to 
a form suitable for these particular patterns and to 
the use of engineering units. 

The potentiometric model cannot be used to obtain 
the constant C appearing in the flow-capacity 
equation of Table 1 because it is not practical 
to use an electrode which is an exact analog of 
the transformed well. However, by using a circular 
electrode scaled to the radius7,,, mentioned pre- 
viously, it is possible to compare measured potential 
drops with calculated values — using, for example, 
Eq. C-14 (Appendix C) modified appropriately to 
apply to the electrical model system. 
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COMPARISON WITH OTHER 
METHODS OF ANALYSIS 


The direct and staggered line drives are special 
cases of the skewed line drive. Areal sweep effi- 
ciencies for these cases have been reported in the 
literature.6,7 In Fig. 5, values computed from the 
analysis developed in this paper are compared with 
values given in Refs. 6 and 7. Fora spacing ratio 
(d/@) = 1.2, the lower limit specified in this paper, 
values for the staggered line drive agree within 1 
areal-sweep-efficiency per cent. The agreement 
for the direct line drive is even better, and both 
cases improve as (d/@) increases. Since these 
cases represent the extremes of skewness, the 
analytical method of this paper should be suitable 
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FIG, 3 —- FLOOD-FRONT ADVANCE AND AREAL SWEEP EFFICIENCY FOR PATTERN A; CALCULATED E 
= 92.2 PER CENT, EXPERIMENTAL E = 91.2 PER CENT. 
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FIG, 4 — FLOOD-FRONT ADVANCE AND AREAL SWEEP EFFICIENCY FOR PATTERN B; CALCULATED E 
= 91.5 PER CENT, EXPERIMENTAL E = 91.4 PER CENT. 
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FIG. 5 — FLOOD-FRONT ADVANCE AND AREAL SWEEP EFFICIENCY FOR PATTERN C; CALCULATED E 
= 87.2 PER CENT, EXPERIMENTAL E = 86.3 PER CENT. 


CONCLUSIONS 


The following conclusions have been made. 

1. This method of analysis provides a simple 
technique for evaluating the areal sweep efficiency 
at breakthrough and the flow capacity associated 
with flooding a formation with anisotropic permea- 
bility, using a developed, skewed line drive. The 
skewed line drive must be considered because a 
direct or staggered line drive may yield, upon trans- 
formation, a skewed line drive. 

2. In anisotropic formations, rows of wells should 
be oriented along the major axis at permeability as 
nearly as possible to obtain larger sweep. In fields 
which are already developed, this can lead to a 
skewed line drive. 


NOMENCLATURE 


C = constant 
E = areal sweep efficiency 
I, = injection-well pattern element 


K = ratio associated with major axis-minor 
axis permeabilities 


P,, production-well pattern element 


Uv 
) 
" 


= stream function cosine 

= stream function co-ordinate 

= functional notation 

= well spacing within rows 

= major axis of ellipse 

= minor axis of ellipse 

= spacing between rows of wells 
= skewness of a pattern 

focal distance of an ellipse 

= thickness 

= permeability 

= general distance co-ordinate 
= general distance co-ordinate 
= index 

= pressure 

= flow rate, reservoir conditions 
= radius 


5 
Xx 
rd 
a 
b 
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f= 
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= travel time 





= velocity 
= distance co-ordinate (transformed system) 


distance co-ordinate (transformed system) 
stream function 


distance co-ordinate (original system) 
distance co-ordinate (original system) 
angle 


tt 


v 
x 
y 

¥ 
n 
¢ 
6 
mM 


viscosity 
a) 
SUBSCRIPTS 


b = breakthrough 
iw = injection well 


porosity 


m = index 
pw = production well 
Ss = streamline index 
= well 
x,y = X, y Co-ordinate system 


> 
uN 
1" 


n, € co-ordinate system; axes oriented 
along major and minor axes of permea- 
bility, respectively 


SUPERSCRIPTS 


~ = transformed co-ordinate system * 





*Note: (7), g) transforms to (x, ¥); @, ¥) rotates to (x, y). 
Quantities unchanged upon rotation from (%, J) to (x, y) continue 
to carry superscript ~. 
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SPACING RATIO, (4/4) 
FIG. 6 — AREAL-SWEEP-EFFICIENCY COMPARISONS. 
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APPENDIX A 


THE POTENTIAL, STREAM-FUNCTION 
AND TRAVEL-TIME EQUATIONS 


The equation presented in the main body of the 
the paper for the pressure distribution in the strip 
of reservoir —0 <x <+%, — (d/2) <-< (d/2) is 


_ = eer). amt). A-1 
p <2 in [cost 22) cos ( 5 (A-1) 


where 
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and q is flow rate at reservoir conditions (negative 
for injection). 
The stream function is given by 


ery 2x 
cos (FF) —| 


cosh( 
Y=- (4,) *arccos 


cosh eee — cos ore 
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as cos(42h¥), 
q 
it follows immediately that 


__ c2e(252) -[1/eomn EX) 
bk cos ee) / = an 1 


The time of travel* along a streamline s from the 
point (x,y) = (0,0) to the point x, y is given by 


» + (A-5) 
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where the velocity component v, or vy is to be 
evaluated along s. 
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To evaluate v, along a particular streamline s, 
Eq. A-5 may be used to give 


h ad —Ccos =n a as ) 
Ree. 











cos 
ecco: 


Substituting in Eqs. A-7b and A-6b gives 
_ (27x 
> sie 
ae (a/2 7) ( 6 ) 


s =(q/4ahe) ; cos (224)-s 


Q 


dx 





(A-9) 





*The variable transformations used make x, y and @ dimen 
sionless. Hence, formally y% and vy have the dimensions (L2/t), 
and ¢ has dimensions (t/L2). This creates no problem in the 
analysis; the conversion to proper units will be shown later. 
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Similar treatment of Eq. A-Gc results in 
ery 
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may also be further manipulated to give 
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s TQ 


The negative sign preceding the equations of Eq. 
A-10 is proper when q is an injection rate, i.e., 
when g has a negative value. If the derivation has 
been carried out for a row of production wells, the 
direction of integration would have been from (x, y) 
toward (0, 0), with a positive sign preceding the 
final Eq. A-10. 

It should be noted that the equations are appropriate 
only in the region 


|x] < (0/2), |y| s (6/2). 
APPENDIX B 


THE AREAL-SWEEP-EFFICIENCY EQUATION 


The generalized line drive is characterized by 
three distance parameters: @, the spacing within 
rows between like wells; d, the spacing between 
rows of unlike wells; and é, the displacement or 
skewness in the x direction of unlike wells. Since 
the choice of positive x direction is arbitrary, e 
may always be taken as positive, ranging from zero 
for the direct line drive to a maximum value (4/2) 
for the staggered line drive. 

The breakthrough time for the streamline which 
passes through the point (X, d/2) will be the sum of 
times in Regions /; and P; (see Fig. 1). Eq. A-10 
may be used to find the time in each region, provid- 
ing consideration is given in Region P; to (ft) the 
displacement @, and (2) the change in sign of the 
rate q. 

Using the absolute value of q, 
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z) (28%), cosh (55) — 
Nxt, “\ 9g ; cos (2X) 
a 














~2 
~ h 
(ty)p, = 2 Jin 


(B-2) 














(B-3) 


To find the point (X, d/2) which yields the initial 
or minimum breakthrough time, Eq. B-3 is differ- 
entiated with respect to X. 


dt ter) (B-4) 
Qa 


“<——* Z (X) sin( 


where Z(X) is positive and finite for 
[E-(6/2)]<x<(6/2), 


the range of X common to Regions |; and P,. If 
the derivative is set equal to zero, 


rex) = 0, 2», etc. . (B-Sa) 
x = (F) E58) EA or oo) 


The only value falling in the X range common to 
Regions |; and P, is 


(B-Sc) 


By differentiating Eq. B-4 a second time, it can be 
shown that X = (€/2) gives a minimum time of travel 
between the injector and producer. 

A similar treatment of the travel-time equation for 
streamlines lying in Regions !, and Py will show 
a second minimum at X = (e — a)/2. Travel time, 
however, is longer for this second minimum; thus, 
the actual minimum travel time occurs along the 
streamline passing through the point (x, y) = (e/2, 
d/2). This is the time of initial breakthrough. In 
effect, each injection well displaces fluid toward 
two different production wells and has a different 
time of breakthrough to each producer. When e = 
(G@/2), the special case of the staggered line drive 
holds, and the times of breakthrough to the two 
producers are equal. 

Substituting Eq. B-5Sc in Eq. B-3 and using the 
even nature of the cosine function, the time of 
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initial breakthrough is given by 


- (498"h) E in/2ia/ al) 





}, "sa eas - 
sti Site cos (7/2) (é/4) sasiias 
The areal sweep efficiency is given by 
__ . oe & @ @ @ -¢ . . . . (B-7a) 
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To complete the working equations, introduce 


K = Ska/ke ie ee ee se ee 


Then, 
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(e/a)’ is not necessarily within the range of 0 to 
0.5. To find (e/a), apply the following where n is 
an integer. If 


n s (&/3)' s (n+ 05), 
then 


oe BT ee eee 


If 
(n-0.5) < (€/3)'< n, 


then 


(6€/G) =n -(6/G)'. . 2. . we ee (Bel Ib) 


Breakthrough time in the original reservoir system 
will be related to t, by 


“) Ste eter 


which may also be written as 
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Areal sweep efficiency will be unchanged from that 
found for the transformed system. 
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APPENDIX C 


THE FLOW-CAPACITY EQUATION 


The flow-capacity equation is developed by con- 
sidering pressure drops in three regions: (1) from 
the transformed-injector wellbore surface (an ellipse) 
to the fictitious circular surface at 7,,, (2) from this 
surface through the main body of the transformed 
reservoir to the isobaric surface located midway 
bet ween the row of injectors and the row of producers, 
and (3) from this surface to the transformed production 
well (also an ellipse). The term 7,,, introduced only 
as an aid in the formal treatment, disappears in 
the final form of the equation. 


THE WELLBORE REGION 


A well centered at the point (0,0) in the original 
(n, ) system of co-ordinates may be described by 
the equation 

2 2 2 
i i ee 


In the x, y co-ordinate system, 
Pe 2 
i. * My / Ky bie. ace 


Yw =O, /ke ’ 


Substituting, the equation describing the transformed 
well is obtained. 
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Xw Yw - 
(7A * Wee tt CD 





This is the equation for an ellipse having semi- 


major and semi-minor axes b,, and C,, given by 


es ee 
2 eo 


The axes of the ellipse are oriented in the ¢ and 7 
directions. 

The transformed system has an equivalent isotropic 
permeability k = VEn . Pressure and flow equa- 
tions are available for isotropic permeability systems 
in which an elliptical source or sink is located in 
an infinite reservoir. These may be applied in the 
immediate vicinity of the transformed well without 
serious error. Muskat® points out that the isobars 
for such a system are confocal ellipses. The focal 
distance is determined by 
































or, since one isobar must coincide with the trans- 
formed well surface, 


As confocal ellipses successively further removed 
from the wellbore are considered, it is clear that 
b and @ approach the same value because / is fixed. 
Stated differently, the isobars become circular in 
form as distance from the wellbore increases. This 
is typical of the pressure distribution due to any 
localized source; at distances sufficiently removed 
from it, the particular geometric form of the source 
no longer influences the pressure distribution. 

The flow equation for the system described is 


; 2 whk(p,,-Py) 


B i= 
Dw + Cw 
where m refers to any general isobar (an ellipse) 


outside the transformed wellbore. At sufficiently 
great distances, 





Eqs. C-4, C-6 and C-7 may be combined to give 
2mhk(p, - 

q = ii 
- In [tArw ete] 





where (r,)esf is an effective well radius given by 
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THE MAIN RESERVOIR REGION 
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It previously has been shown that, essentially, 
an isobar exists along the line y =(d/2), for (d/a) 
> 1.2; i.e., the pressure is independent of x. This 
pressure p(g/2)is obtained from Eq. A-1 by neglect- 
ing the cosine term 


. fe). (22.) ee 
Pavey Gant ae ee 
For (@/@) > 1.2, it follows that 
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For sufficiently small values of 7,,, compared to 
a, the hyperbolic cosine and cosine terms appearing 
in Eq. A-1 may be expanded in series. Only the 
first term of the expansion is significant, so 
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Thus, 
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THE FINAL FLOW-CAPACITY EQUATION 


Eqs. C-8 and C-14 may be combined to eliminate 
Dim and Tm 


(PazziPa) (a) 


(28) +2 in G4- In 2"}} 


This is the pressure drop between y = (d/2) and a 
production well, considering q to be positive. The 
total pressure drop between an injection well and 
production well will be twice this. 
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Using Eq. C-16, 
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It may be noted also that 
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Study of Gas Reservoirs Subject to Water Drive on 


Electronic Differential Analyzer 
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ABSTRACT 


The behavior of gas-storage reservoirs subject to 
water drive is investigated through analog simulation 
on an electronic differential analyzer. The simula- 
tion technique developed on an LM-10 computer per- 
mits the prediction of reservoir volume or pressure 
resulting from the movement of water in the sur- 
rounding aquifer. 

The method developed on the analog computer 
consists of setting up an appropriate transfer-function 
circuit and feeding the arbitrary time-varying bound- 
ary conditions as an input signal. The input may be 
specified as gas reservoir pressure, pore volume or 
the water flux. 

Several cases studied include an isolated gas 
reservoir on a limited aquifer, interference among 
three reservoirs adjacent to a common aquifer and 
the growth of gas-storage volume on an aquifer. 

It is concluded that the method developed on an 
electronic differential analyzer provides an excel- 
lent technique to simulate and investigate the be- 
havior of gas reservoirs subject to water drive. The 
agreement between the reservoir performance as 
predicted from the simulation technique and as 
measured from actual field data is found to be better 
than the range usually encountered in predicting 
water-drive bebavior. 


INTRODUCTION 


It is generally known that some gas-storage res- 
ervoirs are located on top of blanket sands of large 
extent, saturated with brine called aquifers. 

Because the volume of the body of water associ- 
ated with aquifers is usually very large and water 
is compressible, the cyclic pressure variations 
encountered in normal storage service inevitably 
cause unsteady, compressible flow conditions in 
the adjacent aquifers. 

The solution to radial diffusivity equation for a 
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limited aquifer for constant terminal conditions has 
been known! since the early 1930’s. The solution 
for the constant terminal conditions for an infinite 
aquifer was published in 1949.2 

The problem of handling the unsteady flow of 
water through a porous medium with time-varying 
boundary conditions has been studied by many in- 
vestigators through the use of graphical,3 digital4 
computing techniques or on an R-C type of electrical 
analog called a reservoir analyzer.> The digital 
computer performs a series of arithmetic calculations 
to superimpose the solutions of constant terminal 
conditions, while the R-C type of analyzer permits 
the simulation of the continuous medium by means 
of “‘lumped-cell’’ electrical resistance and capaci- 
tance elements based on finite-difference approxi- 
mations. The graphical method is based on the 
convolution integral solution of the partial differ- 
ential equation. 

The purpose of the present study was to simulate 
the behavior of an aquifer on a general-purpose, 
electronic differential analyzer to obtain the influx 
of water or the change in pressure at the reservoir 
continuously and quickly. The use of the electronic 
differential analyzer was based on well known solu- 
tions for the constant terminal conditions already 
available in the literature.2--8 


THEORETICAL DEVELOPMENT 


MODEL, 


The gas reservoir (R) is described in relation to 
an aquifer in Fig. 1. The permeability and porosity 
of an aquifer sand, the viscosity of the brine, and 
the compressibility of the brine and the formation 
are assumed to be constant. It further is assumed 
that the blanket sand is not thick enough to give 
appreciable vertical pressure distribution. The model 
in consideration is circular with radial geometry. 
The aquifer radius can be of limited or infinite ex- 
tent. The gas-water interface is assumed to be fully 
segregated, and the displacement of the gas by 





lReferences given at end of paper. 
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water or water by gas is considered to be pistonlike. 


DIFFERENTIAL EQUATION 


For the model described, the flow of the brine 
through the porous medium can be described by the 
following differential equation. 


o?P, 1 OP _ chp OP (eves te 
@r2 or Or k 00’ 


or in dimensionless form, 


2 P 
A. Spt Ae 


or}, Tp Orp tp 








where 
th . 0 ; 
cours 
and ‘ 
'p = = 
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BOUNDARY CONDITIONS 

The two basic cases encountered in most reser- 
voir-engineering studies are constant-terminal - 
pressure and constant-terminal-rate cases. If the 
pressure is specified at the reservoir boundary, the 
problem may be referred to as terminal-pressure case. 
On the other hand, if the rate is specified one would 
have the terminal-rate case.In gas storage,a variable- 
pressure case often applies; for this case, the pres- 
sure at the aquifer boundary is specified as a func- 
tion of time. 

For the pressure case, the mathematical state- 
ments of the boundary conditions are 


P[1, ty] = Gltp], 
P[oo, ty] = q, 
Plr,,, O] = 0. 


With these boundary conditions, the Laplace trans- 
form of the solution to Eq. 2 is given in the following. 


s+ K, [v3] 


3/2 ‘ge 
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QO. pls]= als] 
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where Q;p(s) is the Laplace transform of the cumu- 
lative water-flux function, and g(s) is the Laplace 
transform of the pressure function G(tp). It can be 
shown that Eq. 3 may be presented in the following 
form 


t 
BD 
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where Q,p is the constant-pressure-case cumulative- 
flux function. 
Eq. 3 can be rewritten to take the form 


Q pls] 
els] 


Eq. 5 becomes a useful equation on an electronic 
differential analyzer and is called a simulation 
equation for terminal-pressure case. Once the Q,){tp) 
function is found, the cumulative water flux or the 
reservoir volume change from the initial reservoir 
pore volume can be calculated by the following 
equation. 


0 
nr kp (9P\ ap 
[2nr 7] “¢] (=) 
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= Sef[sle- +++ ++ 2 ee + (3) 
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Instead of specifying the reservoir pressure, one 
may specify the rate or flux at the reservoir bound- 
ary, with the following boundary conditions. 


oP 
or, [i ty =-E[tp} 


Pleo, t | a O 


i] 


P(r,,0] O. 


The Laplace transform of the pressure-change func- 
tion at the distance of rp can be shown to be 





s¢K, [r,Vs] 
P [rps] = ids : 


s3/2.K [Vs] 


els] {s-/i[rp,s]}- - +--+ (7) 


The terminal-rate case corresponding to Eqs. 4 and 
5 of the terminal-pressure case is given as 


g 
D 
Ra ligiahe f  Pipltp-rar . (8) 
oO T 


where Pp is the constant-rate-case pressure-change 
function, and 
Ap tr ps] 
els] 


The actual pressure change is then 


Pl 7,07) 5): Pen lfnstps- + - « (10) 


As shown in Eqs. 6 and 10, the problem of finding 
the pressure’ change or the reservoir pore-volume 
change narrows down to finding Q,p and P,p. 


s-filrp,s] - eee a) 


The terms /, (s) and {y(p. s) are the Laplace 
transforms of the constant terminal-pressure-case 
cumulative-flux function and the constant terminal- 
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FIG. 1 — RADIAL FLOW MODEL. 


rate-case pressure-change function, respectively. 
Eqs. 5 and 9 are used as the working equations of 
the analog simulation. This analog simulation is 
equivalent to the computation of the superposition 
or convolution integral expressed in Eq. 4 or Eq. 8. 

Unfortunately, the rate data are not available, 
but in most any case the volume-variation data are 
available. From the definition of rate, 


e[s] = i} 1%" Yo! = s-LV,- V, Its] -0. 
dty 


Substituting this equation into Eq. 9, one can obtain 


Pip [rp s] 


- At ae s?filr pSJe ee ee -(11) 
[v, - V, Ils] 
Eq. 11 is the analog simulation equation for the 
volume case. 


ANALOG SIMULATION 


The essential feature of analog simulation on an 
electronic differential analyzer is to obtain the de- 
sired reservoir performance as an output by feeding 
in the boundary condition as an input signal. The 
known information may be the arbitrary reservoir 
pressure schedule, the reservoir'pore-volume varia- 
tion, or the water flux. A block diagram representing 
the sequence of analog simulation is shown on Fig. 
2. 

The term /,(s) in Eq. 5 represents the Laplace 
transform of cumulative-flux funttion of the constant- 
terminal-pressure case or the unit step response for 
the variable-pressure case. Therefore, s/; (s) can be 
approximated from the~solution of the constant- 
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FIG. 3 — INPUT-OUTPUT RELATIONSHIPS FOR THE 
DIFFERENT RESERVOIR-TERMINAL CONDITIONS. 


terminal-pressure case. Since the initial value of 
the cumulative-flux function is zero, from the well 
known property of the Laplace transform rule of a 
derivative, s/,;(s) can be expressed better in the 
following form. 


L IS 4 = s/,[s] -0. 


Fp) 





This relationship can be used to approximate s/,(s) 
from 42s 
dtp 

known, the Eq. 5 can be simulated on an electronic 
differential analyzer. Eq. 5 may be interpreted in 
terms of system transfer function since the transfer 
function is the expression that establishes the re- 
lationship of the input and output variables of a 
system as an operator equation. When constant 
pressure is fed in, the output becomes the constant- 
pressure-case solution. Once the constant-pressure 
response check is made for a linear system, in 
general, the arbitrary pressure response can be 
obtained with the same circuit by simply feeding in 
the proper pressure variation. 10, 11 

The terminal-rate case or volume case can be 
handled by the identical technique on the basis of 
Eq. 9 or Eq. 11 instead of Eq. 5. 

The three different cases of the reservoir terminal 
conditions are represented in Fig. 3. 

When a constant-terminal-case solution is avail- 
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FIG. 2 — OVER-ALL BLOCK DIAGRAM OF THE ANALOG SIMULATION ARRANGEMENT, 
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able, one can find the proper transfer function which 
can be set up on an electronic differential analyzer. 
The desired solution can then be obtained as an 
output, which can be recorded as a function of time. 


EQUIPMENT 


For the studies presented in this paper, the Ster- 
ling model LM 10 electronic differential analyzer 
was used. The LM 10 electronic differential analyzer 
is a small computer equipped with 10 high-gain DC 
amplifiers and 10 ten-turn precision helipots. The 
function G(tp) is generated manually with a 10-turn 
precision helipot. The potentiometer setting is 
changed as a function of time to give the desired 
voltage by means of a vernier dial. The output is 
recorded by using a two-channel Brush recorder. 
Fig. 4 shows a photograph of the actual equipment 
set-up. 


APPLICATIONS 


To demonstrate the application of the theory, a 
few typical problems are chosen. 


PRESSURE CASE FOR AN 
AQUIFER OF LIMITED EXTENT 


As a first illustration and proof of the method, an 
aquifer of rp = 10 is chosen. The rp = 10 is a ratio 
of the size of the aquifer to that of the gas bubble. 
For this aquifer size, three different checks or 
comparisons are made against theoretical or digitally 
computed values. These checks are (1) the unit- 
step-response check with that of van Everdingen and 
Hurst2 table, (2) the application to the actual Field 
A to check with the digital computer solution and 
(3) the comparison of the response for sinusoidal- 
pressure-variation case with that of a mathematical 
solution. 

The first step toward analog simulation is to 
obtain the transfer function based on the previous 
knowledge of unit impulse response. The unit step 
response for rp = 10 is already expressed in an 
equation form, rather than in a tabular form. This 
unit-step response solution for the pressure case 
has been reported for a limited aquifer as follows.? 









FIG. 4 — ELECTRONIC DIFFERENTIAL ANALYZERS 
ARRANGED FOR PRESSURE CASE WITH A FAULT IN 
THE AQUIFER. 
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Since the series converges very rapidly for rp = 10, 
the first term alone gives a good approximation, and 
the flux function for the constant-terminal-pressure 
case for the aquifer size of rp = 10 can be written 
as 


o% = 49.5 = 47 .3e70-0122tp 


= 49.5 _ 47.3¢e70-0122a0 is ax ee (12) 


where a is a constant relating the actual time in 
months to the dimensionless time thy; OF, 


t k0” 3 





sill 2 
COUT. 


and a has a unit of 1/month. It is necessary to 
change the dimensionless time to the actual time 
because the machine is running in the actual time, 
but in a much-condensed scale. For example, one 
month of actual reservoir time is one second of 
machine running time if the machine-time constant 
K is one. In Eq, 12, it is important to note that Q/p 
is dimensionless while Q,p in Eq. 6 has a dimension 
of pressure in pounds per square inch. Upon taking 
the Laplace transform with respect to 6 and multi- 
plying by s, Eq. 12 gives 


47.3s 


sfyls] = anaarre <2" 
S + 0.0122a 


By rearranging, this equation can be written as 
follows. 


sf{s] = 2.2 |1+——213_ 


s 


1 +——_— 
0.0122a 


- (13) 


Based on Eq. 13, an LM 10 circuit diagram for rp = 
10 to get cumulative-flux function is obtained as 
shown in Fig. 5. 

By feeding in a constant voltage at Potentiometer 
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FIG. 5 — LM 10 (ELECTRONIC DIFFERENTIAL ANA- 
LYZER) CIRCUIT DIAGRAM FOR rp = 10 TO GET CU- 
MULATIVE-FLUX FUNCTION, 
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No. 8, which is equivalent to having a constant 
reservoir pressure, the response is obtained. The 
response is the cumulative-flux function at constant- 
pressure case, and the analog result is compared 
with theoretical values? in Fig. 6. The comparison 
in Fig. 6 shows good agreement, and the unit-step- 
response check is essentially a calibration opera- 
tion. It is possible to calibrate the response by re- 
adjusting the potentiometer settings. 

Since the circuit has been established by the unit- 
step-response check, the same circuit can be used 
for any pressure variations by simply generating the 
specified pressure variation with Potentiometer No. 
8. Fig. 7 shows the generated pressure input and 
the function to be generated for the actual Field A. 
For this pressure variation, the response is obtained 
and compared with the results of digital computa- 
tions* in Fig. 8. The reservoir volume is plotted 
against time. 

To secure further proof of the developed method, 
the case of sinusoidal input pressure was studied. 
Table 1 shows the comparison of the analog results 
with theoretical solution,!2 for annual sinusoidal 
pressure variation. 


G = 50 sin (=4- 50 sinl0.0697 tp] psi. 
6 


With the aid of the circuit given in Fig. 5, the 
various effects of practical interest are studied. 
These include the effect of permeability on water- 
return ability for a limited aquifer, the effect of 
aquifer size on water-return ability and the effect 
of phase shift in the pressure cycle on water-return 
ability. Figs. 9, 10 and 11 show the results of this 
study. 
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TABLE 1 — COMPARISON OF ANALOG RESULTS WITH 
THEORETICAL CUMULATIVE-FLUX FUNCTION FOR SINUS- 
OIDAL PRESSURE VARIATION (tp = 7.52 4). 





6 Qip (psi), Quip (psi) 
(month) theoretical Analog Result 
6 633 (efflux) 617 
12 —267 (influx) —275 
24 —356.5 —348 
36 —386.5 —374 
48 —396.2 —396 
60 —399.4 —400 
120 —401 o 





INTERFERENCE AMONG TWO OR MORE 
FIELDS SITUATED ON A COMMON AQUIFER 

The first study of interference among the oil 
fields was made on an electrical reservoir analyzer 
to study the East Texas field.13 But the first math- 
ematical interference study was made by Mortada.® 
The application of Mortada’s method to gas-storage 
problems was made by Coats.14 

The basic idea of this theory stems from the fact 
that the original partial differential equation de- 
scribing the liquid flow through a porous medium is 
linear and homogeneous. 

The pressure drop at a field can be broken down 
into the contributions attributed by the flux of each 
individual neighboring field, as well as by its own 
flux, because the water flux can only be caused at 
the expense of the pressure change. The particular 
example studied in this work was concerned with 
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FIG. 8 — COMPARISON OF ANALOG RESULTS WITH 
DIGITAL RESULTS FOR FIELD A (rp = 10). 
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FIG. 9 — EFFECT OF PERMEAPRILITY ON WATER- 
RETURN ABILITY FOR LIMITED AQUIFER (rp = 10). 














three fields breathing on a common aquifer. 

Fields C and D were discovered in Dec., 1944, 
while Field A was discovered in Sept., 1941. The 
approximate geometry is sketched in Fig. 12. The 
radii of the three fields were approximated from the 
contour maps to be 1.43 to 1.55 miles. The analog- 
predicted volume variation of Field A considering 
the interference of Fields C and D is shown in Fig. 
12 along with actual data. Coats 14 studied Field A 


as an isolated field sitting on a limited aquifer, 
The circuits used for the study are shown in Figs. 
13 and 14, 


MOVING BOUNDARY PROBLEMS 


The problems treated thus far were fixed boundary 
problems. The gas-field radius has been considered 
to be stationary or constant. However, in aquifer 
storage operations, the gas-bubble radius grows from 
zero to a large value during a certain period of time 
before a well developed and stabilized gas bubble 
is established. 

The rigorous analytical solutions applicable to 
aquifer gas-storage problems are just about hopeless 
because of analytical difficulties introduced by the 
moving boundaries. An approximate solution has 
been derived for a growing disk model of a finite 


aquifer thickness. 
ve) 
e dr 
4cdn3/2 pb? [tp - 
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2e + 2e +...f/@r.. (14) 


Eq. 14 has the form of Eq. 8, for which analog sim- 
ulation method is developed. If one defines the unit- 
volume step response for the pressure-change func- 
tion as P,, Eq. 14 can be written 


P 1 p> - Pott.-da 
P = , t,—-T\ aT. 
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FIG. 10 — EFFECT OF AQUIFER SIZE ON WATER-RETURN ABILITY. 
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The unit-volume step response P,7; is computed and method is used. 


is given in Table 2. A circuit to get pressure-change function, using 
In obtaining the pressure-case solution from the the volume as input for the growing disk model of a 

rate- or volume-case solution, the developed analog finite aquifer thickness, is given in Fig. 15. Some 
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FIG, 12—-RESERVOIR-PORE-VOLUME VARIATION FOR FIELD A CONSIDERING THE INTERFERENCE OF FIELDS 
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TABLE 2 — UNIT-VOLUME STEP RESPONSE TO GET 
PRESSURE-CHANGE FUNCTION FOR GROWING DISK OF 
FINITE AQUIFER THICKNESS. 








Pip o? 
tp (Eq. Ra) (Analog Check) 
0. 1 3 1 6 > 
0.5 3.59 i 
1.0 1.773 - 
5.0 0.354 - 
10 0.1772 a 
50 0.03545 - 
100 0.01773 _ 
200 0.00877 _ 
500 0.003545 0.0037 
800 0.00218 - 
1000 0.001815 0,00185 
1300 0.00138 - 
2000 0.00104 0.00105 
3000 0.000782 0.0008 
4000 0.000647 ~ 





corrections have been made to the potentiometer 
settings of P, and P,. With Ps = 0.1952 and P, = 
0.1338, the unit step response for the constant- 
volume case is checked by the analog method and 
is added to Table 2 for comparison. With the same 
circuit of Fig. 15, the constant-pressure-case solu- 
tion is obtained. In doing this, the input (volume 
variation) is manually controlled in such a manner 
that the output pressure-change function would be a 
constant at all times. The constant-pressure-case 
solution thus obtained for the growing disk model 
is shown in Fig. 16. The Q;p in Fig. 16 is the 
quantity defined earlier as Qjp, and upon multiply- 
ing Q;p by a constant, the volume of the gas bubble 
can be obtained. The linearity of the curve at the 
large time portion particularly must be noted. 


Suppose one can assume that the average value 
of the gas-bubble radius times the pressure gradient 
at the gas-water interface is a constant, or 


K=5 correction on P, 
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FIG. 13 — LM 10 CIRCUIT DIAGRAM TO 
GET CUMULATIVE-FLUX FUNCTION FOR 
INFINITE AQUIFER. 





ae 





oP’ i 
"s = Cy dors 
dé /z-0 


Again assuming linear-system response and the gas- 
cap shape to be parabolic, the following relation- 
ship is derived. 


Vi = A Se hese oes 


where 5S, is the time integral of pressure. This 
equation demands that a plot of Qjp°’" vs tp is a 
straight line (Fig. 16). This equation also says that 
a plot of Vz, vs Sy gives a straight line on log-log 
paper. Fig. 17 shows this plot for Field B. The 
slope of the curve in this case should be 4/3, 
theoretically, or 1.333; instead, the slope is found 
to be 1.37 from the actual data correlation. It gives 
a remarkable correlation, however, when the small 
time data are excluded. 

Once the correlation curve or the performance 
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FIG. 14 — LM 10 CIRCUIT DIAGRAM TO 
GET PRESSURE-CHANGE FUNCTION AT 
'pD = 15. 
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FIG.15—LM 10 CIRCUIT DIAGRAM TO GET PRESSURE- 

CHANGE FUNCTION USING VOLUME AS INPUT FOR 

GROWING DISK MODEL OF FINITE AQUIFER THICK- 
NESS, 
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FIG. 16 — CONSTANT-PRESSURE-CASE SOLUTION TO GROWING DISK MODEL OF FINITE AQUIFER THICKNESS. 
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curve is established, as shown in Fig. 17, the 
future predictions can easily be made based on this 
curve. By knowing the gas injection-production 
schedule from the 30th week, the predictions of the 
pore-volume variation and the pressure variation are 
made using Fig. 17; and the predicted results are 
plotted along with the actual data in Figs. 18 and 
19, respectively. The volume variation is predicted 
within 4 per cent error. This small error is magni- 
fied for the pressure prediction, giving as high an 
error as 50 psi when the reservoir pressure is over 


1,000 psia. 


OTHER APPLICATIONS 


The use of an electronic differential analyzer is 
not limited to these examples, but can be extended 
to other reservoir-engineering problems. With slight 
modifications, the present method can be employed 
for oilfield problems. Recently, Hurst!® has pub- 
lished a new chart to handle the oilfield pressure 
decline for constant oil-production-rate case. The 
arbitrary oil production can be handled using the 
new Hurst chart !© by setting up a transfer function 
on an electronic differential analyzer. 

A careful study reveals that there is an alternate 
method which does not necessitate the new chart. 
The usual constant-rate-case solution can be used 
to handle the arbitrary oil-production-rate case in 
simulating an aquifer-oilfield system on an elec- 
tronic differential analyzer. 

For an undersaturated oil field, the aquifer water 
flux can be incorporated with the material-balance 
equation in the following manner. 


Cumulative Water Flux = (ANp —- BP) 


where A and B are constants and Np is the cumula- 
tive oil production. The simulation of this equation 
to obtain the pressure decline from the oil-production 
data is shown in Fig. 20. In Fig. 20, the circuit 
that goes in the block representation is the same 
as that of volume case. The extension of this idea 
can be applied to saturated oil reservoirs, consider- 
ing the gas liberated from solution. 

The method can also be applied to oilfield inter- 
ference studies. A general circuit for oilfield inter- 
ference studies is given in Fig. 21. 
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FIG. 17 — AQUIFER GAS-STORAGE PERFORMANCE 
CURVE FOR FIELD B. 











CONCLUSIONS 


The following general conclusions can be drawn 
from the, study. 

1. The developed, analog simulation method of 
an aquifer on an electronic differential analyzer is 
an excellext method to study the transient behavior 
of an underground gas-storage reservoir. It operates 
continuously, quickly and economically. The accu- 
racy is well within the range usually required for 
typical gas-storage operations. 

2. Because of the difficulties involved in the 
mathematical formulation of the actual problems, 
and the large number of reservoir variables, there 
is no reason why one should resort to an a priori 
model concept based on a series of assumptions. 
The model study can only be a guide to select a 
reasonable unit-step-response curve of a reservoir 
system. 

3. The developed analog simulation method can 
be applied to any other physical linear system to 
analyze its transient behavior. 


NOMENCLATURE 
a = constant relating tp = 20 or A . 
1/month cours 


a 
It 


n = roots of characteristic equation? 


c = sum of aquifer water and formation com- 
pressibility, vol/vol/atm 
c’= sum of aquifer water and formation com- 


pressibility, vol/vol/psi 


e(s) = Laplace transform of E 
E = rate function 


f(s) = Laplace transform of cumulative - flux 
function for constant-pressure case 


{;(s) = Laplace transform of pressure-change 
function of constant-rate case 
g = Laplace transform of G 
G = pressure function 
hb = aquifer water-sand thickness, ft 


Jo. Jy = Bessel functions of the first kind, of 
zero and first order, respectively 


k = permeability, darcies 
K = computer time constant 
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FIG. 18 — COMPARISON BETWEEN FIELD B DATA 
AND PREDICTED PORE VOLUME. 
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FIG. 19 —- COMPARISON BETWEEN OBSERVED AND 
PREDICTED PRESSURES FOR FIELD B. 


Ko, Ky = Bessel function of second kind of zero 
and first order, respectively 


L = Laplace transform 
P = pressure 
P, = pressure at the gas field 
P,‘p = dimensionless pressure-change function 
for constant-terminal-rate case tabu- 
lated in the literature?» ©-8 
P;p = pressure-change function with arbitrary 


terminal rate, cu ft/month 
tD = gtowing-disk-model, unit-volume step re- 
sponse for pressure change, dimen- 
sionless 
(=) = vertical pressure gradient at z = 0 for 
OZ /._6 unit-step pressure, atm/ft 


z= 

QOfp = dimensionless cumulative-flux function 
for constant-pressure case, tabulated 
in the literature?>® 7 

Q yp = cumulative-flux function with arbitrary 


pressure, psi 

5 74 . . . 

Q;p = growing-disk-model, unit-pressure step 
response for cumulative flux, dimen- 
less 

r = radius of aquifer, cm 


, 


r° = radius of aquifer, ft 


r, = gas-field radius 
Tp = 7/r,, dimensionless radius 

s = Laplace transform variable 
S, = time integral of pressure, psi/week 
tp = dimensionless time (see also a) 
Ve = reservoir pore volume, MMcf 
V, = initial reservoir pore volume, MMcf 
W = cumulative water flux, cu ft 

Z = Cartesian co-ordinate 

pi = viscosity of water, cp 

6 = time, month 

6’ = time, second 


= a time variable 


~ 
| 


fractional void space or porosity 


ACKNOWL EDGMENT 


The authors wish to acknowledge the financial 
assistance and the support provided by the Michigan 
Gas Assn. through their fellowship program at the 
U. of Michigan. Valuable assistance, advice and 


SOCIETY OF PETROLEUM ENGINEERS JOURNAL 








al 


1e 
id 








Cumulative 
oil production 


A eee. 











; Volume case 
s* f'(s) transfer 


be function 
‘= l 




















Reservoir pressure 
oo 


FIG. 20 — UNDERSATURATED OILFIELD PRESSURE 
DECLINE FROM OIL-PRODUCTION DATA 


constructive criticism were generously provided by 
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The author applies the theory of hydrodynamic 
stability to fluid flow in porous media and concludes 
that the displacement of miscible liquids is stable 
if the stability coefficient, as defined by Eq. 19.2 
or Eq. 27.0, is positive. The conclusion is invalid, 
however, since it is based on a faulty argument. 

To demonstrate this, the reader is referred to 
Eqs. 4.1 through 4.5 of the paper which describe, 
in differential form, the first-order disturbances c, 
u, v, w and p. In these same equations, Cy, Cz, 21, 
g3 and p are functions of the stable concentration 
c. This concentration ¢ is known and depends on 
x, z and ¢t (the author assumes cy = 0). It follows 
that Fes Cus 21» &3 and 7 are also functions of x, z 
and ¢. With this in mind, we turn to Eqs. 8.1 through 
8.5 which were derived from Eqs. 4.1 through 4.5 
by substituting Eq. 6.0. The significance of this 
substitution is that the author thereby limits the 
solutions of Eqs. 4.1 through 4.5 to equations of 
the form of Eq. 6.0. Eqs. 8.1 through 8.5 represent 
an infinite linear system with five unknowns: c, u, 
v, w and p. The coefficients Cc, Cz, 21, g3 and p 
in this system are known at each point and at each 
instant; the coefficients containing derivatives of / 
are unknown. If this infinite linear system in C, u, 
v, w and p is to have a nontrivial solution, then it 
is indeed necessary (as the author states) that the 
determinant of the coefficient matrix of Eqs. 8.1 
through 8.5, i.e., |A|, be equal to zero. 

This condition, however, is not sufficient. A non- 
trivial solution exists if, and only if, the rank of 
the coefficient matrix of the entire infinite system 
is smaller than five. In other words, all minors of 
order five must vanish and only if all the coefficients 
in Eqs. 8.1 through 8.5 are constant is this equiva- 
lent to the condition |A| =0. Constant coefficients 
means that g,, g3 and p/, for instance, are constant; 
and, since they depend on Cc it follows that ¢ must 
be constant. We conclude that the condition |A| = 
0, on which the author bases his calculations, is 
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necessary and sufficient for a nontrivial solution 
only if the concentration c is constant. In other 
words, the author does not investigate the stability 
of a miscible displacement but, rather, only the 
stability of the flow of a homogeneous liquid. 

Using the initial condition and the fact that the 
coefficients are constant, we find 


fk =@f,=B, fz=nm E=F=G=0. 
Substituting this in Eq. 9.0, the latter becomes 
if, = -a*D, a +n mn”) D3. 
This may be integrated 


f=ax+By+nmnz + ila? D, + (B? + 2 2?)D; It, 


and the stability coefficient 

Yr = a7D, + (B? + n?n*)D3. 
In this last equation, D, and D3 are positive and 
the conclusion is that the motion of a homogeneous 
liquid is stable for all disturbances of the form 

n= nexpi(ax+By+nmz). 
It is perhaps well to point out that at ¢ = 0 these 
disturbances do not satisfy the initial condition of 


Eq. 11.0, which states 


w = 0 for z = 0 and z = 1. 


Therefore, they are unsuitable for investigating the 
stability in a porous medium bounded in the z direc- 
tion by parallel planes (although it would be simple 
to satisfy Eq. 11.0 for w by placing the boundaries 
at z = + 4%, which would result in nonzero values of 


c, at the boundaries and violate the other part of 
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Eq. 11.0). 

The general conclusion is that the author has 
not succeeded in establishing a stability criterion 
for miscible displacement. The reason for this is 
to be found in Eqs. 6.0 and 7.0. By assuming that 
Eqs. 4.1 through 4.5 could be satisfied by Eq. 6.0, 
he automatically excluded the possibility that the 
concentration C varies, and by Eq. 7.0 he excluded 


AUTHOR’S REPLY 


Outmans has suggested, quite correctly, that the 
author’s paper! contains an error. To see what is 
involved requires that one consider both the purpose 
of the study and the methods that are used. 

The comments by Outmans suggest that the Fourier 
series terms (Eq. 6) are not compatible with the 
required boundary conditions. Any real disturbance 
will satisfy flow equations and proper boundary 
conditions. The existence of such a disturbance is 
postulated. Single Fourier terms represent only part 
of the disturbance and, singly, need not satisfy 
such conditions. As a very simple example, consider 
(cy sin mz + c3 sin 3 7 z). For cy =— 3c3, the con- 
dition that the z-direction concentration gradient 
must vanish at the boundaries is satisfied by the 
sum, although not satisfied by either single term. 
Corresponding sine velocity terms vanish identically 
at the boundaries as required. It should have been 
made clear in the original paper that the Fourier 
representation of any real disturbance may include 
such relationships. 

We conclude that any real concentration perturba- 
tion can be described by a Fourier series in the 
space variables x, y and z, with Fourier coefficients 
that will depend only on time. It is important to note 
that the first-order perturbation equations (Eqs. 4.2 
through 4.5) contain time only as a parameter. Thus, 
for a particular Fourier mode for which the spatial 
variation of concentration is known, the set of four 
equations (Eqs. 4.2 through 4.5) can be solved for 
the quantities u, v, w and p. The solutions contain 
a parametric time-dependence. In terms of the pres- 
sure function corresponding to the Fourier mode, 
the set reduces to a linear, second-order, elliptic, 
partial differential equation in canonical form. The 
remaining Eq. 4.1 can be solved to give the time- 
dependence for the Fourier coefficient once u and 
w are determined. Eq. 4.1 for the particular mode 
can be made linear and first-order in the time- 
dependent Fourier coefficient, with the space vari- 
ables as parameters. This approach can be applied 
to inhomogeneous flow systems and is the subject 
of current research effort. 

Most unfortunately, the nature of the pressure 
equation is such that simple solutions for interesting 
systems do not result from the approach just de- 
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flow between parallel z-boundaries from his inves- 
tigation. 

If in any future attempt to find stability criteria 
for miscible displacement in porous media it should 
be necessary to solve Eqs. 4.1 through 4.5, then 
more general solutions than Eq. 6.0 will have to be 
postulated and initial disturbances should be intro- 
duced compatible with the boundary conditions. 


TO H. D. OUTMANS 


scribed. Thus, recourse has been made to other 
methods. Outmans has pointed out a flaw in the ar- 
guments used in one such attempt. However, useful 
results can be achieved in a way which satisfies 
the conditions of rank noted by Outmans. 

Assuming only that the stable flow is well- 
behaved, the coefficients in the first-order pertur- 
bation Eqs. 4.1 through 4.5 of the original paper 
can be expanded as a Taylor series about the posi- 
tion (X%5, Yo» Zo), and about a time defined as ¢t = 0. 
As an example, the coefficient g, in Eq. 4.3 becomes 


Og, og 
81(%,y,2,t) = 27 +|—| (x -— x9) + —} (yY-Yo) 
Ox dy 
oO 
St e~mpelet ts 
Oz - fe] ca 


The superscripts o indicate evaluation of the func- 
tions at (x,,Y%o»Z,0). System properties with a well- 
behaved flow may change rapidly, but never in- 
finitely rapid. Thus, within a sufficiently small 
space region surrounding (x%),¥9,Z 9) and for small- 
enough times ¢, one can neglect all but the leading 
constant term in the series. 

Using standard methods of solution, results are 
obtained that are valid within this limited region 
and which should be compared with those from the 
original paper. One finds that the quantity y, in Eq. 
19.1 disappears and that only the first, time- 
independent term of yp remains in Eq. 19.2. Thus, 
the final result applied to the solution of engineering 
problems,! Eq. 27, is unchanged. This treatment 
amounts to a complete linearization of the stability 
problem and can be extended to inhomogeneous 
systems. The method leads to a prediction of the 
initial growth of a disturbance in the space region 
immediately surrounding a point of interest (Xo, Yo, 
Zo). Thus, while the results are limited, which is 
typical of perturbation theory, they are sufficient 
to satisfy the purpose for which the original de- 
velopment was intended. * & 
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tion of detectors, column types, 
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buys nothing he does not need — 
he chooses options which help his 
work best. For your widest choice, 
call Barber-Colman. 
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packed column work 
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electronic and column 
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23C | 300°C Wy heated filament or thermistor detector 
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™ — Yj  ]Abé/é,yy pl Single or dual packed or 
10 400°C Yj Yj p~ dual capillary columns 

20 400°C GF G7 linear splitter 




















*Convertible to dual column in field. 
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the SJvecti frites recorder 


preferred for laboratory use 


The recti/riter recorder has become the accepted 
laboratory recorder—is preferred for the exacting tasks 
of laboratory applications. The portable recti/riter is 
the only galvanometric rectilinear recorder designed 
specifically as a bench-top instrument with all routine 
controls and adjustments located up front for extra 
convenience. The “writing desk” chart carriage permits 
operators to make the extensive notations usually associ- 
ated with laboratory use while the instrument is recording. 

Ruggedized die-cast construction results in an in- 
strument that can “take it”—yet removal of the one- 
piece dust cover makes every working part completely 
accessible and removable without further disassembly. 
Every recti/riter carries a one-year full service warranty. 

There is a recti/riter to fit your particular require- 
ments—single and dual channel, portable and flush- 
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APPARATUS DIVISION 
PLANTS IN HOUSTON 


AND DALLAS, TEXAS 





mounting models . . . each available in the widest 
selection of standard ranges in the industry. 


Two-Cycle Pen Response 
d-c Milliampere Ranges 
a-c Ampere Ranges 

d-c Ampere Range 


¥2 ma to 100 ma 
0.25 amp to 25 amp 
100 mv for use with 
standard shunts 
Expanded Scale a-c Voltage Ranges 
80-130 V, 160-260 V, 320-520 V 
a-c and d-c Voltage Ranges 10 V to 1000 V 
Frequency Ranges 50, 60, 400 cps 


Five-Cycle Pen Response 


d-c Milliampere Ranges...... ..2.5 ma to 125 ma 


Special options and accessories further expand 
the versatility of recti/riter recorders. Write now for 
complete information on this accepted laboratory 
recorder line. 
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